## UTA MOHRING

Performance Analysis and Capacity Planning of Multi-stage Stochastic Order Fulfilment Systems with Levelled Order Release and Order Deadlines

## BAND 97

Wissenschaftliche Berichte des Instituts für Fördertechnik und Logistiksysteme des Karlsruher Instituts für Technologie (KIT)

Uta Mohring

Performance Analysis and Capacity Planning of Multi-stage Stochastic Order Fulfilment Systems with Levelled Order Release and Order Deadlines

#### WISSENSCHAFTLICHE BERICHTE

Institut für Fördertechnik und Logistiksysteme am Karlsruher Institut für Technologie (KIT)

BAND 97

Performance Analysis and Capacity Planning of Multi-stage Stochastic Order Fulfilment Systems with Levelled Order Release and Order Deadlines

by Uta Mohring

Karlsruher Institut für Technologie Institut für Fördertechnik und Logistiksysteme

Performance Analysis and Capacity Planning of Multi-stage Stochastic Order Fulfilment Systems with Levelled Order Release and Order Deadlines

Zur Erlangung des akademischen Grades einer Doktor-Ingenieurin von der KIT-Fakultät für Maschinenbau des Karlsruher Instituts für Technologie (KIT) genehmigte Dissertation

von Uta Mohring, M.Sc.

Tag der mündlichen Prüfung: 15. Februar 2022 Hauptreferent: Prof. Dr.-Ing. Kai Furmans Korreferenten: Prof. Dr.-Ing. Gisela Lanza, Prof. Dr. Raik Stolletz

**Impressum**

Karlsruher Institut für Technologie (KIT) KIT Scientific Publishing Straße am Forum 2 D-76131 Karlsruhe

KIT Scientific Publishing is a registered trademark of Karlsruhe Institute of Technology. Reprint using the book cover is not allowed.

www.ksp.kit.edu

*This document – excluding parts marked otherwise, the cover, pictures and graphs – is licensed under a Creative Commons Attribution 4.0 International License (CC BY 4.0): https://creativecommons.org/licenses/by/4.0/deed.en*

*The cover page is licensed under a Creative Commons Attribution-No Derivatives 4.0 International License (CC BY-ND 4.0): https://creativecommons.org/licenses/by-nd/4.0/deed.en*

Print on Demand 2022 – Gedruckt auf FSC-zertifiziertem Papier

ISSN 0171-2772 ISBN 978-3-7315-1211-0 DOI 10.5445/KSP/1000148098

# **Danksagung**

Die vorliegende Arbeit entstand während meiner Tätigkeit als wissenschaftliche Mitarbeiterin am Institut für Fördertechnik und Logistiksysteme (IFL) des Karlsruher Instituts für Technologie. Ich möchte mich an dieser Stelle bei allen Personen bedanken, die zum Gelingen dieser Arbeit beigetragen haben.

Zunächst danke ich Prof. Dr.-Ing. Kai Furmans, Leiter des IFL, für die Übernahme des Hauptreferats. Er hat mir viele Freiheiten gewährt und mir stets großes Vertrauen entgegengebracht. Dies hat – gepaart mit den kreativen Diskussionen auf Augenhöhe – zum erfolgreichen Abschluss dieser Arbeit beigetragen. Prof. Dr.-Ing. Gisela Lanza und Prof. Dr. Raik Stolletz danke ich für die Übernahme des Korreferats. Univ.-Prof. Dr.-Ing. Dr. h. c. Albert Albers danke ich für die Übernahme des Prüfungsvorsitzes meiner Promotionsprüfung.

Allen aktiven und ehemaligen Kollegen am IFL danke ich für die angenehme und inspirierende Arbeitsatmosphäre und die gegenseitige Unterstützung. Dies hat nicht nur das Erstellen dieser Arbeit erleichtert, sondern auch meinen Lebensabschnitt am IFL zu einer prägenden und unvergesslichen Zeit gemacht. Insbesondere danke ich Dr.-Ing. Katharina Fleischer-Dörr für ihren Zuspruch und Rückhalt in der Anfangsphase meiner Promotion sowie die kritische Korrektur meiner Arbeit. Bei Christoph Jacobi möchte ich mich für die vielen intensiven Diskussionen und sein stets ehrliches und kritisches Feedback bedanken. Er und Julia Fleischmann haben mich stets darin bestärkt und unterstützt, mich für meine Visionen einzusetzen und meinen eigenen Weg zu gehen. Für diesen Rückhalt und alles, was wir in unserer gemeinsamen Zeit am IFL miteinander und voneinander gelernt haben, bin ich sehr dankbar.

Ein ganz besonderer Dank gilt meinen Eltern und meiner Schwester für ihr Vertrauen, ihr Verständnis und ihre grenzenlose Unterstützung auf meinem bisherigen Lebensweg, aber auch für die notwendige Ablenkung. All das hat mir stets viel Kraft und Rückhalt gegeben. Ich widme diese Arbeit meinem Opa (†2019), dessen Ehrgeiz und Disziplin mir ein großes Vorbild sind.

Karlsruhe, Februar 2022 Uta Mohring

# **Kurzfassung**

Kundenorientierte Auftragsbearbeitungsprozesse in Logistik- und Produktionssystemen sind heutzutage mit kontinuierlich steigenden Auftragsvolumina, hohen Kundenanforderungen hinsichtlich kurzfristiger und individueller Lieferfristen und einer stark stochastisch schwankenden Kundennachfrage konfrontiert. Um trotz der volatilen Kundennachfrage eine effiziente Auftragsbearbeitung und die Einhaltung der kundenindividuellen Lieferfristen gewährleisten zu können, muss die Arbeitslast kundenorientierter Auftragsbearbeitungsprozesse auf geeignete Weise geglättet werden. Hopp and Spearman (2004) unterscheiden zur Kompensation von Schwankungen in Produktionssystemen zwischen den Dimensionen Bestand, Zeit und Kapazität. Diese stellen auch einen guten Ausgangspunkt für die Entwicklung von Glättungskonzepten für stochastische, kundenorientierte Bearbeitungsprozesse dar. In dieser Arbeit werden die Potentiale der Dimensionen Zeit und Kapazität in der Strategie der nivellierten Auftragseinlastung zusammengeführt, um die Arbeitslast mehrstufiger, stochastischer Auftragsbearbeitungsprozesse mit kundenindividuellen Fälligkeitsfristen auf taktischer Ebene zeitlich zu glätten. Ziel dieser Arbeit ist (1) die Entwicklung eines Glättungskonzeptes, der so genannten Strategie der nivellierten Auftragseinlastung, (2) die Entwicklung eines zeitdiskreten analytischen Modells zur Leistungsanalyse und (3) die Entwicklung eines Algorithmus zur Kapazitätsplanung unter Gewährleistung bestimmter Leistungsanforderungen für mehrstufige, stochastische Auftragsbearbeitungsprozesse mit nivellierter Auftragseinlastung und kundenindividuellen Fälligkeitsfristen.

Die Strategie der nivellierten Auftragseinlastung zeichnet sich durch die Bereitstellung zeitlich konstanter Kapazitäten für die Auftragsbearbeitung und eine Auftragsbearbeitung gemäß aufsteigender Fälligkeitsfristen aus. Auf diese Weise wird der zeitliche Spielraum jedes Auftrags zwischen dessen Auftragseingang und dessen Fälligkeitsfrist systematisch zur Kompensation der stochastischen Nachfrageschwankungen genutzt. Die verbleibende Variabilität wird in Abhängigkeit der Leistungsanforderungen der Kunden durch die Höhe der bereitgestellten Kapazität kompensiert. Das analytische Modell zur Leistungsanalyse mehrstufiger, stochastischer Auftragsbearbeitungsprozesse mit nivellierter Auftragseinlastung und kundenindividuellen Fälligkeitsfristen bildet die Auftragsbearbeitung als zeitdiskrete Markov-Kette ab und berechnet verschiedene stochastische und deterministische Leistungskenngrößen auf Basis deren asymptotischer Zustandsverteilung. Diese Kenngrößen, wie beispielsweise Durchsatz, Servicegrad, Auslastung, Anzahl Lost Sales sowie Zeitpuffer und Rückstandsdauer eines Auftrags, ermöglichen eine umfassende und exakte Leistungsanalyse von mehrstufigen, stochastischen Auftragsbearbeitungsprozessen mit nivellierter Auftragseinlastung und kundenindividuellen Fälligkeitsfristen. Der Zusammenhang zwischen der bereitgestellten Kapazität und der damit erreichbaren Leistungsfähigkeit kann nicht explizit durch eine mathematische Gleichung beschrieben werden, sondern ist implizit durch das analytische Modell gegeben. Daher ist das Entscheidungsproblem der Kapazitätsplanung unter Gewährleistung bestimmter Leistungsanforderungen ein Blackbox-Optimierungsproblem. Die problemspezifischen Konfigurationen der Blackbox-Optimierungsalgorithmen Mesh Adaptive Direct Search und Surrogate Optimisation Integer ermöglichen eine zielgerichtete Bestimmung des minimalen prozessspezifischen Kapazitätsbedarfs, der zur Gewährleistung der Leistungsanforderungen der Kunden bereitgestellt werden muss. Diese werden anhand einer oder mehrerer Kenngrößen des Auftragsbearbeitungsprozesses spezifiziert.

Numerische Untersuchungen zur Beurteilung der Leistungsfähigkeit der Strategie der nivellierten Auftragseinlastung zeigen, dass in Systemen mit einer Auslastung größer als 0,6 durch den Einsatz der Strategie der nivellierten Auftragseinlastung ein deutlich höherer α- und β-Servicegrad erreicht werden kann als mit First come first serve. Außerdem ist der Kapazitätsbedarf zur Gewährleistung eines bestimmten α-Servicegrads bei Einsatz der Strategie der nivellierten Auftragseinlastung höchstens so hoch wie bei Einsatz von First come first serve.

# **Abstract**

Order fulfilment systems are forced to manage a highly volatile customer demand consisting of low-volume orders effectively and efficiently while simultaneously meeting customer-required, short order deadlines. Hopp and Spearman (2004) provide three buffer types – inventory buffer, time buffer, and capacity buffer – to handle the volatile workload in production systems. These buffer types also provide several potentials for workload balancing in order fulfilment. In this thesis, we combine the potentials of time buffer and capacity buffer to develop and analytically investigate the Strategy of Levelled Order Release to balance workload in multi-stage, stochastic order fulfilment systems with customerrequired order deadlines over time on a tactical level. The contribution of this thesis is three-fold: We develop (1) a workload balancing concept, the socalled Strategy of Levelled Order Release, (2) a discrete-time analytical model for performance analysis, and (3) an algorithm for capacity planning under performance constraints in multi-stage, stochastic order fulfilment systems with levelled order release and customer-required order deadlines.

The Strategy of Levelled Order Release is characterised by (1) a fixed capacity reserved for order processing in each time period, and (2) an order processing according to ascending due dates in each time period. In this way, the time buffer of each order between its time of arrival and its deadline is used to balance the variability of the customer demand. The remaining variability is compensated by using the capacity buffer depending on the performance requirements of the customers. The developed analytical model for performance analysis of multi-stage, stochastic order fulfilment systems with levelled order release and customer-required order deadlines models system behaviour of such order fulfilment systems as a discrete-time Markov chain. It calculates multiple stochastic and deterministic, system- and customer-related performance measures of order fulfilment systems based on the limiting distribution of the Markov chain. These performance measures, such as system throughput, service level, utilisation, number of lost sales, and time buffer and backlog duration of a completed order, enable a comprehensive and exact performance analysis of multi-stage, stochastic order fulfilment systems with levelled order release and customer-required order deadlines. The relationship between the provided capacity and the performance that is achieved with this capacity cannot be specified by a mathematical equation, but it is given by the analytical model. Thus, the decision problem of capacity planning in multi-stage order fulfilment systems with performance requirements is a blackbox optimisation problem. The problem-specific configurations of the blackbox optimisation algorithms Mesh Adaptive Direct Search and Surrogate Optimisation Integer enable a targetoriented determination of the minimum required, process-specific capacity to meet any performance requirement of the customers that is specified based on one or multiple performance measures of the order fulfilment system.

Numerical studies on the performance evaluation of the Strategy of Levelled Order Release show that one can achieve significantly higher values of α- and β-service level when using the Strategy of Levelled Order Release instead of First come first serve in order fulfilment systems with a utilisation higher than 0.6. Furthermore, the total required capacity to guarantee a predefined α-service level in the order fulfilment system when using the Strategy of Levelled Order Release is at most as high as the one when using First come first serve.

# **Contents**







# **1 Introduction**

Intensified competition due to globalisation and e-commerce increases the pressure on order fulfilment systems to operate in an effective and efficient way (Van Gils et al. 2018; Kundu et al. 2020). Additionally, total order volume constantly increases, while the character of orders has changed from a small number of high-volume orders to a large number of low-volume orders (De Koster et al. 2007). At German Amazon warehouses, for example, an average customer order consists of 1.6 items (Boysen et al. 2019), and in production systems, the trend from mass products to customised products increases the relevance of low-volume/high-variety make-to-order production systems (Kundu et al. 2020). Moreover, promising short, individual, and reliable order deadlines has become a crucial competitive requirement in order fulfilment due to increasing customer orientation (Öner-Közen and Minner 2017). Individual order deadlines depend, for example, on the shipping distance and the service type chosen by the customer (Yan et al. 2010), and delivery promises, such as next-day or same-day delivery, are widespread in e-commerce (Yaman et al. 2012; Boysen et al. 2019). Furthermore, customer demand in order fulfilment systems is highly volatile, independent of the considered industry. Figure 1.1 shows an exemplary time series of incoming customer orders in a German company of the automotive aftermarket sector in 2015. The corresponding squared coefficient of variation of order income is 0.3. The high variability of the number of incoming customer orders results in strong fluctuations of system workload over time (Hopp and Spearman 2004; Boysen et al. 2019). In conclusion, order fulfilment systems are forced to manage a highly volatile customer demand consisting of low-volume orders in an effective and efficient way while simultaneously meeting customer-required, short order deadlines.

Figure 1.1: Time series of incoming customer orders of a German company of the automotive aftermarket sector in 2015 (Gaps in the time series correspond to weekends and public holidays).

Investigating production systems, Hopp and Spearman (2004) claim that "all variability [...] will be buffered". They provide three buffer types to handle the volatile workload in production systems:


In make-to-stock production systems, system workload and customer demand are decoupled by the finished-goods-inventory, so that inventory buffer provides several potentials to compensate the variability of customer demand. In contrast, in the main application fields of order fulfilment – make-to-order production systems and warehouses –, system workload is directly driven by the customer demand since the orders are customer-specific and unknown in advance. Apart from some preparatory tasks, for instance, pre-picking, it is not possible to process the orders before they arrive at the order fulfilment system. Thus, in order fulfilment systems, inventory buffer only provides limited potential to handle the variability of customer demand. Despite the short lead time requirements of customer orders, there is some time buffer in order fulfilment since the processing time of an order is still much shorter than its lead time. However, it is impossible to flexibly adapt order deadlines depending on the current workload since they are given by the customers. Capacity buffer, such as a certain proportion of additional permanent workers or temporary hired workers, provides safety capacity to handle temporary peaks of customer demand. However, using capacity buffer faces the following limitations in order fulfilment: First, due to space and equipment constraints, there is an upper limit on the number of workers that can work simultaneously at the same processing step of the order fulfilment system. Second, due to the increasing probability of blocking, a higher number of assigned workers does not necessarily increase the overall processing performance (Furmans et al. 2009). Finally, providing safety capacity is expensive and might become a competitive drawback (Van Gils et al. 2018; Kundu et al. 2020). In conclusion, we state that inventory, time, and capacity buffer provide multiple, but limited potentials to handle volatile customer demand in order fulfilment systems.

In this thesis, we combine the potentials of time buffer and capacity buffer to develop and analytically investigate the Strategy of Levelled Order Release to balance workload in multi-stage, stochastic order fulfilment systems with customer-required order deadlines. The focus is on the two main application fields of order fulfilment: Warehouses and make-to-order production systems. The contribution of this thesis is three-fold and comprises


## **1.1 Problem Description**

Following the contribution of the thesis, we differentiate three major research segments of the analysis of multi-stage, stochastic order fulfilment systems, resulting in three research questions.

Workload Balancing in Order Fulfilment Systems Following the concept of variability buffers of Hopp and Spearman (2004), time buffer and capacity buffer are suitable instruments to balance variability in order fulfilment systems. By combining the potentials of these buffers, we develop a workload balancing concept for order fulfilment systems that balances system workload over time on a tactical level. The concept is inspired by the concept of Heijunka levelling. Hence, our workload balancing approach is called Strategy of Levelled Order Release. The related research question is the following:

How can we balance workload over time in multi-stage, stochastic order fulfilment systems with customer-required order deadlines?

Performance Analysis of Order Fulfilment Systems To compare the Strategy of Levelled Order Release with alternative strategies and to show the benefit of workload balancing in order fulfilment systems, we model system behaviour of multi-stage, stochastic order fulfilment systems with levelled order release and customer-required order deadlines in a suitable analytical model. Based on this analytical model, we derive exact formulas for various stochastic and deterministic, system- and customer-related performance measures of order fulfilment systems. The related research question is the following:

How can we determine the performance of multi-stage, stochastic order fulfilment systems with levelled order release and customerrequired order deadlines?

Capacity Planning in Order Fulfilment Systems The analytical model enables the performance analysis of any order fulfilment system with a given system configuration. However, the focus of operations managers is not on the performance of a given system configuration but on adapting the current system configuration to guarantee promised performance requirements to their customers. In particular, they have to decide on the capacity that is provided at the processing steps of the order fulfilment system to meet the performance requirements of the customers. For this purpose, we provide an algorithm to solve the capacity planning problem while meeting performance requirements in order fulfilment systems based on the developed analytical model. The related research question is the following:

How can we determine the capacity required to meet specific performance requirements in multi-stage, stochastic order fulfilment systems with levelled order release and customerrequired order deadlines?

## **1.2 Scope of the Thesis**

This thesis is divided into thirteen chapters. Figure 1.2 depicts the overall structure of the thesis and assigns the chapters to the research segments of the thesis.

Initially, Chapter 1 provides the motivation and the problem description of this thesis. Chapter 2 gives a literature review on related topics, such as workload balancing, performance analysis, and capacity planning in order fulfilment systems, and states the research gap of this thesis.

The first research segment on workload balancing in order fulfilment systems comprises Chapters 3 and 4: Chapter 3 introduces a general, formal description of order fulfilment systems that forms the basis for all concepts and models developed in the subsequent chapters. In Chapter 4, we develop the Strategy of Levelled Order Release to balance the workload of order fulfilment systems over time based on the fundamental ideas of Heijunka levelling in production systems and the specific characteristics and requirements of order fulfilment.

The second research segment on performance analysis of order fulfilment systems comprises Chapters 5 to 8: Chapter 5 provides methodological fundamentals and problem-specific specifications of the chosen modelling approach. In Chapters 6 and 7, we introduce an exact and a simplified analytical model for performance analysis of multi-stage, stochastic order fulfilment systems with levelled order release and customer-required order deadlines, respectively. In Chapter 8, we compare these analytical models regarding modelling accuracy, accuracy of performance analysis, and memory and computation time requirements.

The third research segment on capacity planning in order fulfilment systems is discussed in Chapter 9. In this chapter, we formulate the decision problem of capacity planning in order fulfilment systems as a mathematical optimisation model and propose two algorithms to solve the capacity planning problem based on the analytical models introduced in the second research segment.

The last part of this thesis incorporates several numerical studies applying the models and algorithms developed in this thesis: In Chapter 10, we evaluate multiple approaches of runtime and memory optimisation to reduce computation time and memory usage of the analytical models. In Chapter 11, we evaluate the solution algorithms proposed for capacity planning regarding solution quality and runtime efficiency. In Chapter 12, we evaluate the performance of the Strategy of Levelled Order Release in order fulfilment systems compared to alternative strategies.

Finally, Chapter 13 summarises the contribution and the main results of this thesis and presents an outlook on future research directions.

1. Introduction

2. Literature Review

#### **Workload Balancing in Order Fulfilment Systems**

3. Order Fulfilment Systems

4. Strategy of Levelled Order Release

#### **Performance Analysis of Order Fulfilment Systems**


7. Simplified Model for Performance Analysis

8. Evaluation of Models for Performance Analysis

#### **Capacity Planning in Order Fulfilment Systems**

9. Formalisation and Solution Algorithms of the Capacity Planning Problem

#### **Numerical Analysis and Application**

10. Runtime and Memory Optimisation of the Markov Chain


#### **Conclusion and Outlook**

13. Conclusion

# **2 Literature Review**

This chapter presents the state of the art of academic literature in the three research segments of this thesis: Workload balancing in order fulfilment systems in Section 2.1, analytical models for performance analysis of order fulfilment systems in Section 2.2, and capacity planning approaches in the context of order fulfilment in Section 2.3. The focus of each section is on the two main application fields of order fulfilment: Make-to-order production systems and warehouses. The state of the art of academic literature results from a keywordbased systematic literature review in the Scopus database. The methodology of the literature review and the keywords are described in Appendix A. The results of the literature review are summarised in Section 2.4 by defining the research gap of this thesis.

## **2.1 Workload Balancing**

Workload balancing is defined as managing the variability of system workload over a specific time horizon (Irastorza and Deane 1974). Academic literature discusses workload balancing in different application fields: Project scheduling, vehicle routing, facility balancing, assembly line balancing, and Heijunka levelling. Table 2.1 provides a classification of workload balancing approaches in these application fields.

In project scheduling, resource capacity is fixed, and there are precedence constraints between the project tasks. The decision problem is to determine the start


Table 2.1: Workload balancing approaches in different application fields (extended representation based on Vanheusden et al. (2020)).

times of the project tasks provided that maximum project duration and precedence constraints are met, and provided that resource utilisation is balanced over time (Rieck et al. 2012). Vehicle routing determines the routes of a given set of vehicles so that customer demand is fulfilled within the required time windows, and vehicles' capacity utilisation is balanced. Matl et al. (2018) provide a review of workload balancing approaches in vehicle routing. In facility balancing, the workload is balanced among multiple facilities and over time by adequately assigning a given set of independent tasks to a given set of facilities (Huang et al. 2006). Assembly line balancing problems assign tasks to stations of the assembly line so that the total workload for manufacturing one product unit is balanced among the stations and provided that precedence constraints between the tasks and the fixed number of available machines are met (Becker and Scholl 2006). Heijunka levelling balances both production volume and product mix of a multi-product production line over time. The decision problem of Heijunka levelling is on determining the buffer size of finished-goods-inventory that is required to guarantee a predefined service level (Matzka et al. 2012).

Workload balancing in order fulfilment differs from the above mentioned balancing approaches to some extent (see Table 2.1): First, resource capacity in order fulfilment is flexible since a high proportion of tasks is still manual (De Koster et al. 2007; Marchet et al. 2015; Peeters and Van Ooijen 2020), and it is possible to quickly hire additional temporary workers in case of shortage (Van Gils et al. 2017). Hence, we initially balance system workload over time. Subsequently, we determine the capacity that is required to meet the performance requirements. Second, the tasks in order fulfilment are time-constrained by clear, individual, and externally given order deadlines (Öner-Közen and Minner 2017).

In the following, we provide the state of the art of workload balancing in the two main application fields of order fulfilment: Make-to-order systems in Section 2.1.1 and warehouses in Section 2.1.2.

### **2.1.1 Workload Balancing in Make-To-Order Systems**

The classification of order review and release (ORR) policies of Bergamaschi et al. (1997) mentions workload balancing as a possible decision criterion of workload control in make-to-order production systems. Balanced workload control focuses on balancing the workload among the workstations of the production system. Optimisation approaches are used to minimise the deviations of the current workload from its target workload level at every workstation (Kundu et al. 2020). This lean-based approach differs from the majority of ORR policies in make-to-order production systems, which predominantly focus on strictly limiting the workload of every workstation (Kundu and Staudacher 2017). Kundu et al. (2020) provide a literature review of ORR strategies incorporating workload balancing.

While balanced workload control in job shops is widely discussed in the literature, there are only few publications on balanced workload control in flow shops: Portioli-Staudacher and Tantardini (2012) introduce a lean-based ORR policy for flow shops that aims at levelling workload along the flow. The proposed optimisation model minimises the unbalancing of workload among workstations and the unbalancing over time. Fernandes et al. (2020) develop an optimisationbased ORR policy that combines characteristics of balancing and limiting order release approaches. The corresponding optimisation model minimises the differences between the released workload and a predefined workload norm at every workstation so that lower and upper bounds of the stations' workload are met. Kundu et al. (2020) introduce a modular design of ORR policies, conduct a simulation study to better understand the interaction effects between the parameters of an ORR policy in flow shops, and develop new ORR policies for flow shops with high processing time variability. The simulation results show that an ORR policy which uses an optimisation approach to balance the workload among the workstations enables higher system performance than ORR policies focusing on limiting the workload at the workstations.

The order fulfilment systems studied in this thesis can be seen as flow shops. However, the above mentioned publications on balanced workload control differ from our workload balancing concept regarding the balancing unit since we study workload balancing over time and not among workstations (see Table 2.2).


Table 2.2: Workload balancing concepts in the context of order fulfilment in make-to-order systems (above), in warehouses (middle), and in this thesis (below) (extended representation based on Vanheusden et al. (2020)).

<sup>1</sup> The concept of variability buffers only applies to workload balancing concepts using the balancing unit "time".

### **2.1.2 Workload Balancing in Warehouses**

In warehousing, workload balancing is applied in different settings (see Table 2.2). Zone picking is a classical approach to reduce the time needed for order processing in order picking systems. For this purpose, the order picking system is subdivided into multiple zones, and each order can be simultaneously processed in the different zones. Jane (2000) and Jane and Laih (2005) provide two approaches for workload balancing among the zones of a zone order picking system in the long term: Zone sizing and zone assignment. Jane (2000) adjusts the zone size in order to balance the workload among the pickers, whereas Jane and Laih (2005) develop a clustering algorithm that assigns items to zones in order to balance the workload among the pickers. Dynamic picking systems, such as bucket brigades, consider workload balancing among workers in manual order picking systems in the short term. In the concept of bucket brigade, multiple pickers work in a flow line, whereby every worker follows one simple rule: Go on processing an order until your successor in the flow line is met, who takes over the order. Then, go back upstream and take over the current order of your predecessor in the flow line. If the pickers are sequenced according to their picking performance from the slowest to the fastest, the workload of the order picking system is balanced, independent of the starting positions of the workers (Bartholdi III and Eisenstein 1996). Academic literature studies different approaches for throughput improvement of bucket brigades and different configurations of order picking systems. Hong (2019) provides a current literature review on bucket brigades. In contrast to the above mentioned workload balancing concepts, this thesis focuses on workload balancing over time (see Table 2.2).

Flexible workforce planning represents a reasonable approach to balance temporary fluctuations in customer demand: Workers of a pool of fixed and temporary workers are flexibly scheduled to balance system workload over time and among workers. This approach expects a precise forecast of the daily workload. Based on forecasted workload and worker productivity, the total number of required workers and their allocation is determined. Van Gils et al. (2017) propose this workload balancing approach for a zone order picking system and present forecasting methods to precisely forecast the daily number of incoming order lines. Flexible workforce planning concentrates on the capacity buffer for workload balancing, whereas the Strategy of Levelled Order Release regards a combination of capacity and time buffer (see Table 2.2).

Vanheusden et al. (2020) and Vanheusden et al. (2021) study workforce balancing as a planning step in daily resource capacity planning of warehouses:


The authors introduce the Operational Workload Balancing Problem (OWBP) to evenly distribute the tasks in an order picking system over the time periods of one day, whereby the number of tasks, their release, and their deadline are given. The objective function minimises the unbalancing of the scheduled workload over all time periods of one day, which is measured by the range (Vanheusden et al. 2020), the maximum, the lexicographic maximum, the variance, or the Gini coefficient (Vanheusden et al. 2021) of the scheduled workload. The OWBP is solved by an iterated local search algorithm. This approach differs from the Strategy of Levelled Order Release regarding decision level and the mathematical modelling of system parameters: Vanheusden et al. (2020) and Vanheusden et al. (2021) analyse a planning horizon of one day, such that customer demand and order deadlines are deterministic, and the workload is balanced on an hourly level. In contrast, we examine a planning period of multiple days. Thus, customer demand and order deadlines are stochastic, and the workload is balanced on a daily or a shift-based level (see Table 2.2).

## **2.2 Performance Analysis of Order Fulfilment Systems**

The literature review on analytical models for performance analysis of order fulfilment systems is limited to the two main application fields of order fulfilment: Make-to-order systems in Section 2.2.1 and warehouses in Section 2.2.2.

## **2.2.1 Performance Analysis of Make-To-Order Systems**

Analytical models, especially queueing models, for performance analysis of stochastic manufacturing systems are widely discussed in academic literature. Dallery and Gershwin (1992), Buzacott and Shanthikumar (1993), Papadopoulos and Heavey (1996), Govil and Fu (1999), and Papadopoulos et al. (2019) provide comprehensive reviews on queueing models in manufacturing. However, analytical models that focus particularly on performance analysis of make-toorder production systems are barely studied in the literature. Haskose et al. (2002) and Haskose et al. (2004) introduce queueing models for performance analysis of workload control production planning in make-to-order systems with limited buffer capacities that provide exact results for small arbitrary make-toorder systems. To analyse make-to-order systems of any size and complexity, the authors provide an approximation approach. In contrast, simulation models are commonly used to analyse make-to-order systems (Haskose et al. 2004).

Furthermore, performance analysis of make-to-order systems is often limited to system-related performance measures, such as utilisation, throughput, and inventory levels, although they are interrelated with customer-related performance measures, such as service level and order tardiness (Altendorfer and Jodlbauer 2011). Operations managers in make-to-order systems are often confronted with a trade-off between system- and customer-related performance measures: While high capacity and inventory levels are reasonable from the customer-related perspective to guarantee high service levels, they are adverse from the system-related perspective due to low system utilisation and high operational costs. Although customer-related performance measures are key to quantify the perceived quality of service, they are often neglected in analytical models for performance analysis of make-to-order systems (Altendorfer and Jodlbauer 2011).

Customer service level measures to which extent the delivery time requirements of the customer orders are met. Thus, knowing the order deadlines is key when calculating customer service level. In the research fields of order acceptance and due date setting in make-to-order systems, order deadlines are modelled as endogenous variables: In the job entry phase of workload control, customer enquiry management determines whether to bid for an order or not and, if so, what due date and price should be (Thürer et al. 2011). Common strategies of customer enquiry management are total acceptance, acceptance based on present and future workload, and due date negotiation (Mezzogori et al. 2021), whereby the focus of academic literature is on the latter two. These two strategies have in common that orders can be rejected, and that order deadlines are determined internally. Slotnick (2011) and Thürer et al. (2019) provide literature reviews of order acceptance and due date setting, respectively. In contrast, order deadlines that result from the shipping distance, the service type chosen by the customer or delivery promises, such as next-day or same-day delivery (Yan et al. 2010; Yaman et al. 2012) are exogenous variables. These customer-required order deadlines are predominantly defined to be deterministic and constant for every order (Altendorfer 2014). Liu and Yuan (2001) specify order deadlines by a deterministic maximum delivery lead time in their queueing model for performance analysis of a simple assembly network where different components are processed at two workstations before being merged at the assembly station. Assuming exponentially distributed interarrival times, processing times, and machine breakdowns, the authors provide exact formulas for the probability distribution of the flow time, the throughput, the expected work-in-process, and the service level, and introduce a mathematical program to maximise system throughput while maintaining the required customer service level. Alternatively, Souza and Ketzenberg (2002) model order deadlines by a deterministic upper bound of the average order lead time when analysing a two-stage make-to-order remanufacturing system in order to determine the long-run production mix that maximises profit subject to a service level constraint.

However, real make-to-order systems face random, customer-required order deadlines due to different products ordered by different customers and the planning uncertainty the customers of a manufacturing company face (Altendorfer and Minner 2011). There are only few publications regarding both stochastic, customer-required order deadlines and customer-related performance measures in make-to-order systems: The heavy traffic analysis of single-stage stochastic systems with the scheduling policy Earliest due date, conducted by Doytchinov et al. (2001), incorporates generally distributed customer-required order deadlines, interarrival times, and processing times. The authors provide approximate formulas of customer service level and the probability distribution of customer lateness. Moreover, the priority sequencing problem in make-to-order systems with respect to minimum tardiness cost, which is modelled as Markovian decision process by Öner-Közen and Minner (2017), incorporates stochastic, customer-required order deadlines as well as the customer-related performance measures service level and expected order tardiness. The model is primarily used to evaluate the impact of priority sequencing decisions on the performance measures. However, exact performance analysis of single-stage, stochastic make-to-order systems can be conducted based on the model when using a simple due date-based priority sequencing policy, such as Earliest due date. Altendorfer and Jodlbauer (2011) introduce an analytical model for performance analysis of single-stage, stochastic make-to-order systems with stochastic, customer-required order deadlines that incorporates multiple customer-related performance measures, such as service level, expected lead time, expected finished-goods-inventory, and expected order tardiness. Assuming exponentially distributed interarrival times, processing times, and order deadlines, the computed performance measures are exact. For systems with arbitrary distributed parameters, the performance analysis is limited to approximations. Altendorfer and Minner (2011) extend this model to two-stage make-to-order systems in order to analyse the simultaneous optimisation of manufacturing capacities and planned lead times with respect to total inventory holding and customer order tardiness cost.

Table 2.3 summarises the relevant analytical models for performance analysis of stochastic make-to-order systems by classifying them regarding modelling approach, modelling of order deadlines, number of stages, performance measures, and solution quality. To the best of our knowledge, the model of Altendorfer and Minner (2011) is the only one for performance analysis of multi-stage make-toorder systems incorporating stochastic, customer-required order deadlines and customer-related performance measures. However, it differs from our model regarding the modelling approach: Altendorfer and Minner (2011) provide a continuous-time model whose solution quality depends on the shape of the probability distributions of the stochastic parameters, whereas our discrete-time model provides exact results for any arbitrary distributed parameters. Furthermore, discrete-time models enable the calculation of complete probability distributions of relevant stochastic performance measures (Schleyer and Furmans 2007).

### **2.2.2 Performance Analysis of Warehouses**

The focus of research in warehousing is on the analysis and optimisation of specific planning and control problems, especially in order picking, such as batching, routing, zoning, sequencing, worker assignment, and storage assignment. There are numerous analytical models analysing the impact of specific control strategies on warehouse performance and determining optimal parameter values and strategy configurations in order to maximise warehouse performance. Rouwenhorst et al. (2000), Gu et al. (2007), De Koster et al. (2007), and Davarzani and Norrman (2015) provide comprehensive reviews on warehouse operations. Despite the existing interdependencies between different processes and decisions in warehouse operations, integrated modelling and analysis of multiple warehouse planning problems as well as a holistic performance analysis of warehouses are barely studied in the literature so far (Davarzani and Norrman


2015; Van Gils et al. 2017). Van Gils et al. (2018) state the potentials of integrated analysis and simultaneous solution approaches of multiple operational planning problems and provide a classification of combined planning problems in warehouse operations. The authors conclude that integrated modelling and analysis focuses on batching and storage assignment, routing and storage assignment, and batching and routing. However, these integrated models regard combinations of multiple specific decision problems. None of them studies a holistic performance analysis of warehouses.

Furthermore, stochastic problems are rarely studied in warehouse literature, although warehouse managers suffer from several stochastic influences and uncertainties resulting from both the outside supply chain and internal processes. Deterministic models can provide good approximations of stable stochastic systems. However, when analysing warehouses with highly fluctuating parameters, deterministic models may strongly deviate from actual system behaviour and lead to wrong conclusions (Gong and De Koster 2011). Gong and De Koster (2011) state the potentials of stochastic models and optimisation approaches and provide an overview of stochastic research in warehouse operations.

The existing stochastic models for performance analysis of warehouses, which especially investigate different configurations of order picking systems, focus on system-related performance measures, such as travel time (Chew and Tang 1999; Le-Duc and De Koster 2007; Pan et al. 2014), throughput time (Tang and Chew 1997; Gong and De Koster 2011; Le-Duc and De Koster 2007; Yu and De Koster 2009; Van Nieuwenhuyse and De Koster 2009; Pan and Wu 2012; Xu et al. 2014), and throughput capacity (Van Der Gaast et al. 2020). The provided models are approximations, and they are restricted to expected values and variances of the performance measures. In contrast, customer-related performance measures, such as service level and order tardiness, are rarely considered despite their importance with respect to customer satisfaction (Van Gils et al. 2018).

Consequently, order deadlines that are essential for quantifying customer-related performance measures are barely incorporated into stochastic warehouse models. In the wave picking models of Çeven and Gue (2017) and MacCarthy et al. (2019), customer orders that arrive before a certain cutoff time are due to the next deadline corresponding to a truck departure in the shipping area of the warehouse. Thus, each wave is specified by a cutoff time and a deadline, and service level is measured as the proportion of orders that is completed until their deadline. Based on this performance measure, called Next Scheduled Deadline, Çeven and Gue (2017) optimise the timing and the number of order waves, and MacCarthy et al. (2019) provide best performance frontiers to determine minimum picking rate as well as optimum timing and number of order waves. Furthermore, the discrete-time model for performance analysis of a one-block warehouse with order batching of Schleyer and Gue (2012) specifies order deadlines by a deterministic maximum throughput time. Based on this model, the authors calculate the probability distribution of the throughput time and the service level and derive the optimal batch size with respect to a service level constraint.


Table 2.4: Analytical models for performance analysis of stochastic systems in warehousing.

Table 2.4 summarises the relevant analytical models for performance analysis of stochastic systems in warehousing. Although the aforementioned models incorporate order deadlines, they are restricted to deterministic order deadlines. Furthermore, they are limited to order picking, that is, single-stage systems. Thus, there is no analytical model in warehouse literature for performance analysis of stochastic, multi-stage order fulfilment processes with stochastic, customer-required order deadlines with respect to customer-related performance measures.

## **2.3 Capacity Planning in Order Fulfilment Systems**

Capacity planning determines the resource requirements to meet a given customer demand over a planning horizon (Chen et al. 2009). Capacity planning problems with customer-required deadlines and due date-related performance constraints are barely studied in the literature. Table 2.5 presents an overview of the related literature that can be subdivided into performance analysis and comparison of capacity planning policies and optimisation of capacity planning. Furthermore, the publications differ regarding the number of stages, the modelling of order deadlines, and the chosen performance metric.

Bertrand and Sridharan (2001) evaluate subcontracting policies in single-stage systems with stochastic, customer-required lead times, Poisson-distributed order arrivals, and negative exponentially distributed processing times. Subcontracting is necessary since the order arrival rate is greater than the service rate. The developed subcontracting policies specify when and which orders should be subcontracted and vary regarding informational needs and complexity. In a numerical study, these policies are compared regarding the percentage of tardy orders, the expected tardiness of an order, capacity utilisation, and throughput time. Mincsovics and Dellaert (2009) propose workload-dependent capacity planning policies in a single-stage system with exponentially distributed interarrival and service times and a fixed lead time. The number of unprocessed orders


Table 2.5: Capacity planning approaches with due date-related performance constraints: Capacity planning policies (above) and optimisation of capacity planning (below).

is constrained by a constant upper bound, and all arriving orders exceeding this bound are refused. The workload-dependent capacity planning policies are defined by two switching points. A switching point is a specific workload value, for which an order arrival triggers a switch in capacity. The authors present one exact and one approximate analytical approach to evaluate the due date-related performance of the policies, measured by the expected number of backorders and the expected number of early orders. A policy with constant provided capacity is used as a benchmark in the numerical study. Altendorfer et al. (2014) investigate periodic capacity planning policies regarding potential improvements of service level and the expected tardiness of an order in multi-stage systems with stochastic customer demand, processing times, and customer-required deadlines. The policies differ regarding the amount of available information on the stochastic system parameters and the approach used to derive the amount of provided capacity from the amount of demanded capacity. A policy with constant provided capacity is used as a benchmark for these policies in a simulation study.

Altendorfer and Minner (2011) and Altendorfer and Minner (2012) investigate capacity planning in multi-stage systems with exponentially distributed interarrival times, processing times, and customer-required deadlines. Altendorfer and Minner (2011) examine the planned lead time release policy that releases orders to the next stage whenever the remaining time to due date is shorter than the sum of planned lead times of the remaining stages. In contrast, Altendorfer and Minner (2012) examine the work-ahead window release policy that releases an order when its remaining time to due date is shorter than the work-ahead window. Using a continuous-time queueing model, Altendorfer and Minner (2011) model the system by a series of M|M|1 queueing systems, whereby the provided capacity per stage is modelled by the stage-dependent processing rate. Altendorfer and Minner (2012) model the system by a series of M|M|m queueing systems and the provided capacity per stage by the processing rate and the number of servers per stage. Based on the queueing model, the authors derive analytic formulas for the expected values of work-in-process per stage, finished-goodsinventory, and the number of backorders. The capacity planning problem is formulated as an unconstrained optimisation problem that minimises the expected costs of work-in-process per stage, finished-goods-inventory, backorders, and capacity. Besides the provided processing rate per stage, further decision variables are considered: The planned lead time per stage in Altendorfer and Minner (2011) and the number of servers per stage and the work-ahead window in Altendorfer and Minner (2012). Ma et al. (2019) study capacity planning with a response time constraint in single-stage systems. Customer demand is Poisson-distributed, processing times are generally distributed, and there is a deterministic maximum response time per order. The system is modelled as a M|G|1 queueing system, whose processing rate specifies the required capacity. Based on the Pollaczek-Khintchine formula, the expected sojourn time of an order is calculated. The response time constraint is modelled by penalty costs that occur if the expected sojourn time of an order exceeds the maximum response time per order. Thus, the capacity planning problem is to minimise the sum of capacity cost and expected penalty cost. Buyukkaramikli et al. (2013) examine capacity planning with a throughput time constraint in single-stage systems with exponentially distributed interarrival times, processing times, and customer-required lead times. Capacity can be changed periodically, whereby two capacity levels are possible: Permanent capacity level and permanent plus contingent capacity level. The throughput time constraint guarantees that each order is completed before its predefined lead time with a certain probability. The optimisation problem for capacity planning determines the period length, the permanent capacity level, the contingent capacity level, and the switching point at which contingent capacity is deployed, such that the capacity-related cost is minimised, and the throughput time constraint is satisfied. For this, the system is modelled as a queueing system, and the periodic capacity policy is reflected in the change of the processing rate. The capacity planning problem is solved using a search algorithm.

Except Buyukkaramikli et al. (2013), the literature on capacity planning with customer-required deadlines and due date-related performance constraints is limited to average-based performance metrics, such as the expected number of backorders, the expected sojourn time, or the expected tardiness (see Table 2.5). In contrast, in this thesis, the performance constraint of capacity planning is specified by the probability distribution of a performance measure, such as system throughput or backlog duration, or the service level of the order fulfilment system.

## **2.4 Chapter Conclusion**

The results of the literature review state the research gap of this thesis. First, there is no workload balancing concept in the literature to balance the workload of order fulfilment systems over time on a tactical level by using a combination of time buffer and capacity buffer (see Table 2.2). Second, the discrete-time analytical model provided in this thesis is the first analytical model for performance analysis of workload balancing in order fulfilment systems that ensures a real-life representation of order fulfilment systems (see Tables 2.3, 2.4). It regards customer-required order deadlines and the interdependencies between the processing steps, and it models customer demand, customer-required deadlines, and processing performance as stochastic parameters. Furthermore, the discrete-time analytical model enables the exact calculation of complete probability distributions for several system- and customer-related performance measures of the order fulfilment system. Third, capacity planning problems with due date-related performance constraints, which rely on probability distributions of corresponding performance measures and service level, in multi-stage systems with stochastic, customer-required order deadlines have not been studied in literature so far (see Table 2.5).

# **3 Order Fulfilment Systems**

This chapter aims at introducing a general and formal description of order fulfilment systems, as they form the basis for all concepts and models developed in the following. In Section 3.1, we give an overview of the definition of order fulfilment in the literature. Section 3.2 introduces a formal description of order fulfilment systems. Section 3.3 summarises the results of this chapter.

## **3.1 Definition of Order Fulfilment**

The Global Supply Chain Forum identifies order fulfilment as one of the eight key processes of supply chain management, along with customer relationship management, customer service management, demand management, manufacturing flow management, procurement, product development and commercialisation, and returns (Cooper et al. 1997). The process of order fulfilment aims at designing a logistics network that permits a company to meet customer requests while minimising total costs. It involves all activities on defining customer requirements, designing the logistics network, and filling customer orders (Croxton 2003). The activities of order fulfilment are subdivided into


Strategic order fulfilment is on establishing a suitable structure for managing the order fulfilment process effectively and efficiently by incorporating the requirements of manufacturing, logistics, and marketing. Table 3.1 gives an overview


Table 3.1: Sub-processes of strategic order fulfilment (adopted from Croxton (2003)).

of the sub-processes of strategic order fulfilment. In contrast, operational order fulfilment focuses on the execution of the order fulfilment process once it has been established. It determines how customer orders are generated and communicated, recorded, processed, documented, picked, and delivered (Croxton et al. 2001; Croxton 2003). Table 3.2 gives an overview of the sub-processes of operational order fulfilment.


Table 3.2: Sub-processes of operational order fulfilment (adopted from Croxton (2003)).

Lin and Shaw (1998) identify the order fulfilment process as one of the three pillars of a company, along with the product development process and the customer service process. The order fulfilment process "starts with receiving orders from the customers and ends with having the finished goods delivered" (Lin and Shaw 1998). It involves order management, manufacturing, and distribution. Order management is on receiving customer orders and committing order requests. Manufacturing involves production scheduling, material planning, capacity planning, and shop floor control, and distribution focuses on inventory and transportation planning. Following Lin and Shaw (1998), the main objectives of order fulfilment are


Okongwu et al. (2012) extend the definitions of Lin and Shaw (1998) and Croxton et al. (2001) by subdividing the order fulfilment process into the following phases:


Order promising verifies resource availability to commit reliable due dates to the customers. In contrast, order fulfilment focuses on the execution of the order fulfilment process after the due dates are determined. It incorporates operational disruptions, such as machine breakdowns, material shortage, and quality defects (Okongwu et al. 2012).

## **3.2 Formal Description**

In this thesis, the focus is on operational order fulfilment. We concentrate on the phase of order fulfilment of the order fulfilment process since order promising is irrelevant under the assumption of customer-required order deadlines.

We describe order fulfilment from a superior, abstract, and process-oriented perspective (see Figure 3.1): Incoming, unprocessed orders with customer-required

Figure 3.1: Exemplary structure of an order fulfilment system with three order types (green, red, blue) and ten processes (A-J).

order deadlines determine the starting point of the order fulfilment process (see definition of Lin and Shaw (1998)). These orders are processed within the order fulfilment system in a sequence of multiple processing steps. We do not further specify the processing steps of the order fulfilment system in this thesis. We abstract from a subdivision into, for instance, administrative, manufacturing, and logistics tasks or a subdivision according to the sub-processes of operational order fulfilment (see Table 3.2). Instead, we model every processing step as a separate process that is specified by the parameters introduced in Section 3.2.3. Completed orders ready for delivery represent the output of the order fulfilment system (see definition of Lin and Shaw (1998)).

To evaluate the performance of the order fulfilment system regarding the overarching objectives of order fulfilment, we calculate several system- and customerrelated performance measures. System-related performance measures characterise the efficiency of the order fulfilment process. Common examples are the utilisation, the number of unprocessed orders, and the number of lost sales. In contrast, customer-related performance measures specify the characteristics of the completed orders and are indicators to quantify customer satisfaction. Common examples of customer-related performance measures are the service level and the backlog duration of a completed order.

The resources that are required for order processing at each processing step of the order fulfilment system are not explicitly modelled as a separate parameter. Instead, they are integrated into the specification of each processing step by determining the time capacity provided for order processing at each process. This approach is sufficient since we abstract from the technical design and implementation of the processing steps – manual, partly automated, or fully automated – and thus from different resource types, such as workers of different qualification levels, different variants of tools and machines, or any combination of multiple resource types. If necessary, the number of required resources at any processing step can be derived from the amount of required time capacity by considering its technical design and its specific characteristics. However, this is beyond the scope of this thesis.

We abstract from the occurrence of material shortage by assuming that all products ordered by the customers are always available in the required quantity in the order fulfilment system. Thus, aspects of inventory management and product availability are neglected in the concepts and models developed in the following.

In order fulfilment, orders are usually classified into different order types. Orders are subdivided, for instance, regarding their deadlines into standard and express orders. Furthermore, in warehouse operations, a classification of the orders in picker-to-parts and parts-to-picker can be meaningful.

In conclusion, we formally define an order fulfilment system by


i ∈ I to an adjacency matrix, whose rows and columns correspond to the processes p ∈ P of the order fulfilment system:

$$\begin{aligned} v: \mathcal{I} &\to \{0, 1\}^{p\_{max} \times p\_{max}} \\ i &\mapsto \mathbf{V}\_i = \begin{pmatrix} v\_{i,1,1} & \dots & v\_{i,1,|\mathcal{P}|} \\ \vdots & \ddots & \vdots \\ v\_{i,|\mathcal{P}|,1} & \dots & v\_{i,|\mathcal{P}|,|\mathcal{P}|} \end{pmatrix} \end{aligned} \tag{3.1}$$

with

$$v\_{i,p,p'} = \begin{cases} 1 & \text{if process } p' \text{ succeeds process } p \text{ in} \\ & \text{processing sequence of order type } i \quad \forall i \in \mathcal{I}, \forall p, p' \in \mathcal{P}. \\ 0 & \text{otherwise} \end{cases} (3.2)$$

#### **3.2.1 Specification of an Order**

An order represents the reference base in order fulfilment. However, an order can be specified in different ways: In the simplest case, every order corresponds to a single customer order. Furthermore, internal processing orders that result from the customer orders depending on the processing strategy of order fulfilment can represent the reference base in order fulfilment. In order batching, customer orders are consolidated to batches by some criteria, whereby every batch represents one internal processing order. Order batching is widespread in logistics and production systems since it is often more efficient to process a batch of customer orders instead of processing a single customer order. For instance, there are processing steps in production systems at which a certain number of units is processed simultaneously. Thus, it is reasonable to group multiple customer orders to one batch to efficiently use the capacity at these processing steps. Moreover, in order picking, multiple customer orders are pooled to one batch to increase the picking performance (Schleyer and Furmans 2007). In contrast, in order picking systems with zone picking, every customer order is subdivided into several internal processing orders, one per zone, to simultaneously process the customer order in the different zones (Gu et al. 2007; De Koster et al. 2007). Furthermore, in the case of high-volume customer orders, it can be meaningful to subdivide the customer orders into their order lines, so that every order line represents one internal processing order.

The concrete specification of the reference base is not relevant for the concepts and models developed in this thesis, provided that all parameters specifying the order fulfilment system refer to the same reference base. For simplicity reasons, we sum all specifications of the reference base under the expression order and only use this expression in the following.

### **3.2.2 Specification of an Order Type**

Each order type i ∈ I of an order fulfilment system is specified by the parameters number of incoming orders and lead time of an order:

Number of incoming orders The number of incoming orders A<sup>i</sup> specifies the volume of orders of order type i that arrive at the order fulfilment system per time period. It is modelled as a discrete random variable with the finite range A<sup>i</sup> = {ai,min, . . . , ai,max} to ensure a real-life representation of the highly volatile customer demand in order fulfilment systems.

Lead time of an order The lead time of an order E<sup>i</sup> specifies the time buffer of an order of order type i between its time of arrival and its deadline. It is measured in number of time periods. We model the lead time of an order as a discrete random variable to incorporate the real-life characteristics of order deadlines in order fulfilment: Individual and customer-required deadlines. The range of the lead time of an order is given by E<sup>i</sup> = {ei,min, . . . , ei,max}. The minimum possible lead time of an order is zero. In this case, the order is due in the same time period in which it arrives at the system.

### **3.2.3 Specification of a Process**

Each process p ∈ P of an order fulfilment system is specified by the parameters processing performance and capacity:

Processing performance The processing performance Li,p specifies the number orders of order type i that can be completely processed per time unit at process p. Operational disruptions, such as machine and tool breakdowns and rework of orders, have an impact on the processing performance. They are incorporated by modelling the processing performance as a discrete random variable. Its range is given by Li,p = {li,p,min, . . . , li,p,max}.

Capacity The capacity ci,p specifies the amount of time capacity that is provided for order processing of order type i at process p in every time period. It is measured in time units, and it is modelled as a deterministic, integer, and non-negative parameter.

## **3.3 Chapter Conclusion**

In this chapter, we presented the formal description of an order fulfilment system that forms the basis for all concepts and models developed in the following. We define an order fulfilment system by a finite set of order types, a finite set of processes, and a function mapping each order type to its processing sequence in the order fulfilment system. Each order type is characterised by the number of incoming orders and the lead time of an order. Each process is specified by its processing performance and its capacity.

# **4 Strategy of Levelled Order Release**

This chapter aims at developing a levelling concept for order fulfilment systems that balances system workload over time. The concept is inspired by the key ideas of Heijunka levelling in production systems. In Sections 4.1 and 4.2, we present the concept of Heijunka levelling in production systems and discuss the differences between levelling in order fulfilment and levelling in production systems. Based on these findings, we develop a levelling concept for order fulfilment systems, the so-called Strategy of Levelled Order Release, in Section 4.3. By summarising the results of this chapter, Section 4.4 provides an answer to the first research question of this thesis.

## **4.1 Concept of Heijunka Levelling in Production Systems**

Section 4.1 is based on Mohring et al. (2020). Parts of the following text and the figures are taken from that publication without changes.

Heijunka levelling is a simple and widespread concept in production systems to manage the production of several product types on one common production line. It has its origins in the Toyota Production System. The key idea of Heijunka levelling is to convert the volatile customer demand into a regular, recurring, and standardised production schedule to guarantee an even utilisation of the given production capacity. Furthermore, Heijunka levelling aims at supplying downstream processes with a constant flow of small lots of different product types and at generating a constant demand of parts from upstream processes. Consequently, the need of excess capacity or stocks to handle peaks as well as the bullwip effect are reduced or even eliminated. Heijunka levelling smooths both the volume and the product mix of the production system (Matzka et al. 2012; Liker 2004, p.115f.).

Heijunka levelling is applied to one selected production step of the production line, the so-called pacemaker. This pacemaker production step controls the production sequence and the production rhythm of the whole production line: Kanban loops control upstream production steps, and downstream production steps form a continuous flow system. The choice of the pacemaker depends on several characteristics of the production line, such as product portfolio, production capacity, and customer requirements. The pacemaker production step is usually close to the customer. However, the further upstream the pacemaker the better (Furmans and Veit 2013; Smalley and Womack 2004, p.30).

The planning procedure of Heijunka levelling refers to one planning period of the production system, for instance, one month or one quarter. The planning period is subdivided into smaller scheduling intervals, such as one week, one day, or one shift. The planning procedure of Heijunka levelling consists of two planning steps: System parametrisation and Operational planning (see Figure 4.1).

Figure 4.1: Planning procedure of Heijunka levelling.

## **4.1.1 System Parametrisation**

System parametrisation occurs at the beginning of every planning period and aims at smoothing both the production volume and the product mix.

Volume smoothing determines the production capacity per product type per scheduling interval that is reserved for producing this product type in each scheduling interval. For this, the total product type-specific customer demand of the planning period is evenly distributed on the scheduling intervals. The reserved product type-specific production capacity corresponds to the average customer demand per scheduling interval of this product type (Matzka et al. 2012; Furmans and Veit 2013).

Product mix smoothing determines the production sequence of the product types within a scheduling interval. Common objectives of production sequence planning are minimising setup times or maximising the regularity of the product mix. These decision problems are discussed in detail by the research area of level scheduling. Dhamala and Kubiak (2005) and Boysen et al. (2009) provide comprehensive reviews of level scheduling. The production sequence is characterised by the so-called EPEI ("each part each interval"). The EPEI corresponds to the time interval in which at least one part of each product type can be produced (Matzka et al. 2012).

Figure 4.2 presents a model of a Heijunka levelled kanban system. The levelling pattern on the Heijunka-board visualises the reserved production capacities and the production sequence per scheduling interval.

## **4.1.2 Operational Planning**

Based on the levelling pattern, operational planning occurs at the beginning of each scheduling interval. The incoming customer orders of the current scheduling interval are completed by taking the required products from the

Figure 4.2: Model of a Heijunka levelled kanban system.

finished-goods-supermarket. The associated kanbans return to the Heijunkaboard. These kanbans are assigned to the reserved production capacities of the corresponding product types on the Heijunka-board according to First come first serve (see Figure 4.2). If the customer demand of a product type exceeds its reserved production capacity in the current scheduling interval, the associated kanbans are stored in an overflow box. They are assigned to future scheduling intervals when the customer demand of this product type is below its reserved production capacity. If the reserved production capacity of a product type exceeds the number of free kanbans of this product type, the excess production capacity remains unused in order to avoid the production of unneeded products (Matzka et al. 2012; Furmans and Veit 2013).

## **4.2 Delimitation from Heijunka Levelling**

The workload balancing concept of Heijunka levelling balances a volatile customer demand over time, as it is the objective of workload balancing in order fulfilment systems. However, due to differences in the system characteristics and the overarching question of the levelling concept, it is impossible to directly apply the concept of Heijunka levelling in order fulfilment systems. Table 4.1 provides an overview of these differences, and they are discussed in detail in the following.


Table 4.1: Comparison of levelling in order fulfilment systems and Heijunka levelling in production systems regarding system characteristics (above) and overarching question (below).

First, production systems and order fulfilment systems differ regarding order lot size. In order fulfilment systems, orders are customer-specific in product type and product volume. The probability that multiple customers place the same order at the same time is negligible. Therefore, order lot size in order fulfilment is usually one. It is possible to have a lot size of one in production systems. However, it is predominantly more efficient to have larger lots in production systems due to the required setups between the production of different product types. Setup processes, such as tool changes, cleaning operations, and programming setups, are usually time-consuming, especially in immature production systems. Consequently, higher lot sizes lead to a decreasing frequency of setup processes and higher productivity in production systems (Dehdari 2014, p.18).

Second, customer demand plays a different role in Heijunka levelling than in the levelling concept for order fulfilment systems since the focus of Heijunka levelling is on make-to-stock production systems. Make-to-stock production systems satisfy the customer demand from the finished-goods-inventory (see Figure 4.2). Thus, system workload and customer demand are decoupled by the finished-goods-inventory in Heijunka levelling. In contrast, in order fulfilment systems, system workload is directly driven by the customer demand since the orders are customer-specific and unknown in advance. Apart from some preparatory tasks, it is not possible to process the orders before they arrive at the order fulfilment system.

Third, there are differences regarding the characteristics of the available capacity. In order fulfilment systems, a high proportion of processing steps is still manual due to the high flexibility of workers and the high investment cost of automated systems (De Koster et al. 2007; Marchet et al. 2015; Peeters and Van Ooijen 2020). Hence, capacity in order fulfilment systems is predominantly determined by the number of assigned workers. Furthermore, capacity is flexible in the short term since it is possible to hire additional temporary workers in case of shortage (Van Gils et al. 2017). In contrast, capacity in production systems is determined by the number of available machines, which is fixed in the short term (Becker and Scholl 2006).

Finally, levelling in order fulfilment systems differs from Heijunka levelling regarding the overarching question of levelling. In both application areas, the objective is to guarantee specific promised performance requirements, such as a service level of 99%, at minimum costs. However, the starting points to achieve this objective are different due to the above mentioned differences in system characteristics: The question of Heijunka levelling is to decide on a suitable buffer size of the finished-goods-inventory to meet the promised performance requirements (Furmans 2005; Lippolt and Furmans 2008; Matzka et al. 2012; Furmans and Veit 2013). In contrast, in order fulfilment systems, the decision is on a suitable amount of provided capacity to meet the promised performance requirements.

## **4.3 Levelling Concept for Order Fulfilment Systems**

Section 4.3 is based on Mohring et al. (2020). Parts of the following text are taken from that publication without changes.

The levelling concept for order fulfilment systems balances the volatile system workload over time and comprises two planning problems: Capacity planning and order dispatching. It is called the Strategy of Levelled Order Release. It is inspired by the key ideas of Heijunka levelling (see Section 4.1) and considers the specific characteristics and requirements of order fulfilment (see Section 4.2).

In order fulfilment systems, the pacemaker processing step always corresponds to the first processing step in the processing sequence of an order type since orders are customer-specific and system workload directly depends on the customer demand (see Section 4.2). The Strategy of Levelled Order Release is applied to the pacemaker. All downstream processing steps form a continuous flow system.

The planning procedure of the Strategy of Levelled Order Release consists of the planning steps of system parametrisation and operational planning, analogous to Heijunka levelling (see Figure 4.1).

## **4.3.1 System Parametrisation**

System parametrisation occurs at the beginning of every planning period and determines the capacity per order type per scheduling interval that is reserved for order processing of this order type in each scheduling interval (smoothing of volume), and the processing sequence of the order types within a scheduling interval (smoothing of product mix).

The reserved capacity per scheduling interval c¯i,p of order type i ∈ I at process p ∈ P is calculated as the ratio of the expected number of incoming orders per scheduling interval E(Ai) of order type ito the expected processing performance per time unit E(Li,p) of order type i at process p:

$$
\vec{c}\_{i,p} = \frac{E(A\_i)}{E(L\_{i,p})}.\tag{4.1}
$$

This approach of capacity planning should be seen as one possible, straightforward approach to determine a good initial capacity level. However, this capacity level is probably insufficient in order fulfilment systems that face high performance requirements, such as service level requirements of 99%. A solution approach to determine the minimum capacity that is required to guarantee specific performance requirements in the order fulfilment system will be presented in Chapter 9.

To determine the processing sequence of the order types, we use the same methodology as in Heijunka levelling (see Section 4.1.1). The resulting levelling pattern is visualised on the Heijunka-board (see Figure 4.3) and forms the starting point of the operational planning.

Figure 4.3: Model of an order fulfilment system with levelled order release.

### **4.3.2 Operational Planning**

Operational planning occurs at the beginning of each scheduling interval and allocates unprocessed orders of the different order types to the corresponding reserved capacities in the levelling pattern. The pool of assignable orders comprises the order income of the current scheduling interval and the unprocessed orders remaining from previous scheduling intervals stored in the overflow box of the Heijunka-board (see Figure 4.3). To determine the processing sequence of orders of a specific order type, we consider their individual due dates as follows:


If the number of unprocessed orders of an order type exceeds its reserved capacity in the current scheduling interval, the order backlog in the overflow box increases by the corresponding number of orders. Otherwise, the remaining capacity is deployed for training, maintenance, and continuous improvement measures.

## **4.4 Chapter Conclusion**

In this chapter, we developed a levelling concept for order fulfilment systems, the so-called Strategy of Levelled Order Release, based on the key ideas of Heijunka levelling in production systems. It is impossible to directly apply the concept of Heijunka levelling in order fulfilment systems due to several differences in the system characteristics between production systems and order fulfilment systems (see Table 4.1). Furthermore, the overarching question of the levelling concept is different: Heijunka levelling decides on a suitable buffer size of the finished-goods-inventory to meet the promised performance requirements, whereas in order fulfilment systems, the decision is on a suitable amount of provided capacity to meet the promised performance requirements.

The key elements of the Strategy of Levelled Order Release are the following:


The Strategy of Levelled Order Release balances volatile workload in multistage, stochastic order fulfilment systems with order deadlines over time by combining the potentials of time buffer and capacity buffer in order fulfilment systems: By processing the orders according to ascending due dates, the time buffer of each order between its time of arrival and its deadline is systematically exploited to balance the variability of the customer demand. The remaining variability of the customer demand is either passed on to the customer, resulting in low service levels, or it is balanced using the capacity buffer. Thus, the Strategy of Levelled Order Release provides an answer to the first research question of the thesis:

### How can we balance workload over time in multi-stage, stochastic order fulfilment systems with customer-required order deadlines?

The extent to which capacity buffer is used to balance the remaining variability depends on the specific performance requirements of the order fulfilment system. According to the Strategy of Levelled Order Release, the reserved capacity per order type per scheduling interval is fixed within one planning period. However, at the beginning of every planning period, capacity can be adapted individually depending on the specific performance requirements of order fulfilment since capacity in order fulfilment systems is flexible in the short term (see Section 4.2). The approach of capacity planning based on expected values presented in Section 4.3.1 should be seen as one possible, straightforward approach to determine a good initial capacity level. This capacity level is probably insufficient in order fulfilment systems that face high performance requirements, such as a service level of 99%. A solution approach to determine the minimum capacity that is required to guarantee specific performance requirements will be presented in Chapter 9.

# **5 Choice and Specification of Modelling Approach**

This chapter aims at selecting and specifying a suitable modelling approach to model and analyse system behaviour and performance of multi-stage, stochastic order fulfilment systems with levelled order release and customer-required order deadlines. We initially discuss the advantages of discrete-time analytical models compared to simulation models and continuous-time models (see Section 5.1) and introduce the theoretical basics of discrete-time Markov chains (see Section 5.2). Subsequently, Section 5.3 provides several specifications for modelling system behaviour of multi-stage order fulfilment systems with levelled order release and customer-required order deadlines as a discrete-time Markov chain. Finally, Section 5.4 summarises the results of this chapter.

## **5.1 Discrete-time Analytical Model**

Analytical models and simulation models are widespread modelling approaches to depict and analyse system behaviour and performance of stochastic systems. Simulation models provide a high degree of freedom and enable modelling at any desired level of detail. However, modelling, validation, and performing experiments using simulation models is time-consuming, and simulation models are usually tailored to the individual characteristics of the considered system. In contrast, analytical models are of higher generality and more time efficient. However, analytical models suffer from a limited level of detail since they are often restricted to simplified representations of the real system. Furthermore, key figures calculated based on simulation results only approximate actual system behaviour, whereas analytical models enable an exact calculation of relevant key figures (Schleyer 2007, p.17-19). Due to the drawbacks of simulation models compared to analytical models, we prefer an analytical model to analyse order fulfilment systems.

Analytical models are differentiated into continuous-time and discrete-time models. Continuous-time models depict system behaviour at any point in time t ∈ R0, whereas in discrete-time models, time is discretised into constant time intervals of length tinc. System behaviour is observed at points in time t that are integer multiples of the discretisation interval tinc (Stewart 1994, p.4f.):

$$t = k \cdot t\_{inc} \qquad k \in \mathbb{N}\_0. \tag{5.1}$$

Discrete-time models provide several advantages regarding accuracy, level of detail, and the description of real processes compared to continuous-time models (Schleyer 2007, p.13-17): First, discrete-time models provide higher accuracy in specifying the input parameters since complete probability distributions are given to specify generally distributed stochastic parameters, whereas in continuous-time models, they are only specified by the first two moments. Second, calculated key figures in continuous-time models are often limited to their expected values and their variances. However, it is essential for practitioners to know specific quantiles, for instance, when calculating the service level. Hence, complete probability distributions of the key figures are required. Discrete-time models enable the exact calculation of complete probability distributions of key figures and thus provide a more detailed analysis of system behaviour than continuous-time models. Third, real processes, such as the number of incoming orders and the processing performance in order fulfilment systems (see Section 3.2), are discrete. Furthermore, probability distributions resulting from real data that is observed in as-is analyses of the real system are discrete and of arbitrary shape. In continuous-time models, these discrete and arbitrary probability distributions are approximated by suitable parametric, theoretical probability distributions, which leads to additional computational effort and an imprecise description of the real processes. In contrast, in discrete-time models, the discrete and arbitrary probability distributions are directly used without any additional computational effort and any loss of precision. Due to these advantages of discrete-time models, we choose a discrete-time analytical model to analyse order fulfilment systems.

When modelling an order fulfilment system as a discrete-time analytical model, one realises that the system state occurring at time (t + 1) only depends on the current system state at time t and the realisations of the stochastic parameters occurring between time t and time (t + 1). This property is called the Markov property. Hence, a discrete-time Markov chain is a suitable model to analyse system behaviour of order fulfilment systems.

## **5.2 Discrete-time Markov Chain**

A discrete-time Markov chain is a stochastic process {X<sup>t</sup> , t ∈ T } with a discrete time parameter set T = {0, 1, . . .} and a discrete state space X = {0, 1, . . .} whose conditional probability function satisfies the Markov property:

$$P(X^{t+1} = i^{t+1} \mid X^0 = i^0, \dots, X^t = i^t) = P(X^{t+1} = i^{t+1} \mid X^t = i^t), \qquad (5.2)$$

whereby X<sup>t</sup> specifies the state of the stochastic process observed at time t. Hence, state X<sup>t</sup> contains all relevant information concerning the history of the process (Stewart 1994, p.4f.).

The conditional probabilities P(X<sup>t</sup>+1 = i <sup>t</sup>+1 | X<sup>t</sup> = i t ) are called single-step transition probabilities or transition probabilities for short. They give the conditional probability of making a transition from state i t at time t to state i <sup>t</sup>+1 at time (t + 1). If the transition probabilities are time-independent, the Markov chain is said to be homogeneous. The transition matrix P summarises the transition probabilities in a (|X | × |X |)-dimensional stochastic matrix, whereby the entry pi,j in the i th row and the j th column corresponds to the probability of making a transition from state i at time t to state j at time (t + 1):

$$p\_{i,j} = P(X^{t+1} = j \mid X^t = i) \qquad i, j \in \mathcal{X}. \tag{5.3}$$

The n-step transition probability

$$p\_{i,j}^{(n)} = P(X^{t+n} = j \mid X^t = i) \qquad i, j \in \mathcal{X},\tag{5.4}$$

is the probability of reaching state j at time (t + n), in n transitions, when starting in state i at time t (Stewart 1994, p.5f.).

#### **5.2.1 Characteristics**

The states of a discrete-time Markov chain and the Markov chain itself have several mathematical characteristics. In the following, we only introduce those characteristics that become relevant in the subsequent chapters of the thesis. The following definitions are based on Privault (2013, p.117-126).

Communicating States A state j ∈ X is said to be accessible from state i ∈ X if it is possible to reach state j with non-zero probability in a finite number of transitions when starting in state i. If states i and j are accessible from each other, respectively, they are said to be communicating states. A subset of the state space X <sup>0</sup> ⊂ X is called a communicating class if every pair of states i, j ∈ X <sup>0</sup> is communicating.

Irreducible Markov Chain A Markov chain whose state space consists of a unique communicating class is said to be an irreducible Markov chain. Otherwise, the Markov chain is called reducible.

Recurrent and Transient States A state i ∈ X is said to be a recurrent state if starting in state i, the Markov chain will return to state i within a finite time with probability of one. Otherwise, in the case of a probability smaller than one, state i is said to be transient. States of the same communicating class are either all recurrent or all transient. A communicating class of recurrent (transient) states is called a recurrent class (transient class).

Aperiodic States The period of a state i ∈ X is the greatest common divisor of n ∈ N with p (n) i,i > 0. State i is said to be aperiodic if it has a period of one. Otherwise, state i is periodic. All states of the same communicating class have the same period. A Markov chain is said to be aperiodic if all states are aperiodic.

### **5.2.2 Relevant Probability Distributions**

When studying the behaviour of a discrete-time Markov chain, the following probability distributions are of particular interest.

State Probability Distribution at Time t The state probability distribution at time t

$$\begin{aligned} \pi^t &= \begin{pmatrix} \pi\_0^t & \pi\_1^t & \dots & \pi\_{\left| \mathcal{X} \right|}^t \end{pmatrix} \\ \pi\_i^t &= P(X^t = i) \qquad \forall i \in \mathcal{X}, \end{aligned} \tag{5.5}$$

specifies the probability of being in any state i ∈ X at time t. Given the initial state probability distribution π 0 , the state probabilities at time t are calculated as follows (Stewart 1994, p.14)

$$
\boldsymbol{\pi}^t = \boldsymbol{\pi}^0 \cdot \mathbf{P}^t. \tag{5.6}
$$

Stationary Distribution State probabilities of a Markov chain are said to be stationary if any transition according to the single-step transition probabilities P have no effect on the state probabilities (Bolch 2006, p.41):

$$
\boldsymbol{\pi} = \boldsymbol{\pi} \cdot \mathbf{P}.\tag{5.7}
$$

In the case of an irreducible Markov chain with finite state space, the stationary distribution π is obtained by solving the following set of linear equations:

$$\pi\_j = \sum\_{i \in \mathcal{X}} \pi\_i \cdot p\_{i,j} \quad \forall j \in \mathcal{X} \tag{5.8}$$

$$\sum\_{i \in \mathcal{X}} \pi\_i = 1.\tag{5.9}$$

In the case of a reducible Markov chain with finite state space, its stationary distribution π consists the stationary distributions of its recurrent classes that are calculated by separately solving the set of linear equations (5.8)-(5.9) for every recurrent class. The stationary probability of every transient state is zero (Cinlar 1975, p.152-156).

Limiting Distribution The limiting distribution specifies the asymptotic behaviour of the Markov chain. Given the initial state probability distribution π 0 of the Markov chain, if the limit

$$\tilde{\pi} = \lim\_{t \to \infty} \pi^t \tag{5.10}$$

exists, then this limit is called the limiting distribution of the Markov chain (Stewart 1994, p.15). For any aperiodic Markov chain, the limiting distribution π˜ exists. In the case of an aperiodic, irreducible Markov chain with finite state space, the limiting distribution π˜ is independent of the initial state probability distribution π 0 and corresponds to the unique stationary distribution π of the Markov chain (Bolch 2006, p.47). In the case of an aperiodic, reducible Markov chain with finite state space, the limiting distribution π˜ depends on the initial state of the Markov chain. If the initial state is recurrent, the limiting distribution π˜ corresponds to the stationary distribution of its recurrent class. Otherwise, if the initial state is transient, the limiting distribution π˜ is calculated as the weighted sum of the stationary distributions of the recurrent classes. The weights correspond to the probabilities that a certain recurrent class is accessible from the initial state. The matrix of limiting distributions Π˜ summarises the limiting distributions π˜ for any initial state, whereby the i th row corresponds to the limiting distribution of initial state i ∈ X (Cinlar 1975, p.152-156).

## **5.3 Modelling Order Fulfilment as a Discrete-time Markov Chain**

We use the general description of an order fulfilment system introduced in Section 3.2 to model system behaviour of multi-stage order fulfilment systems with levelled order release and customer-required order deadlines as a discrete-time Markov chain. Section 5.3.1 introduces several specifications of the object of consideration, and Section 5.3.2 specifies the considered modelling approaches.

### **5.3.1 Specification of the Order Fulfilment System**

We consider the different order types i ∈ I of the order fulfilment system separately since the reserved capacities per time period are order type-specific in the case of levelled order release (see Section 4.3), so that there are no interaction effects between different order types. To analyse an order fulfilment system with multiple order types i ∈ I, we use one separate Markov chain for each order type.

We assume each order type i ∈ I to have a sequential processing sequence and the processes p ∈ P to be sorted according to the processing sequence of the considered order type. Thus, process p = 1 is the first processing step and process pmax is the last processing step of order type i in the order fulfilment system. An order type whose processing sequence contains splits is modelled by dividing this order type into multiple "artificial" order types, each of which has a sequential processing sequence.

We consider orders with failed due dates, but we limit the possible backlog duration of an order by the maximum accepted backlog duration of R time periods. Consequently, orders become lost sales once their backlog duration exceeds the maximum accepted backlog duration. The set of due dates K of an order is limited by the minimum possible due date of an order of (−R) time periods and the maximum possible lead time of an order of emax time periods:

$$\mathcal{K} = \{-R, \dots, e\_{\max}\}.\tag{5.11}$$

### **5.3.2 Specification of the Modelling Approach**

There are different approaches to model the multi-stage character of order processing and the transmission of orders between successive processing steps in a multi-stage order fulfilment system as a discrete-time Markov chain. In this thesis, we provide the following modelling approaches:


The exact modelling approach models each processing step of the multi-stage order fulfilment system as a separate process (see Figure 5.1). In this manner, the exact modelling approach considers the interdependencies between the processing steps by exactly recording the partially processed orders that are transmitted between successive processing steps of the multi-stage order fulfilment system. However, the high level of detail of the exact modelling approach leads to high modelling complexity and high computational effort of the corresponding discrete-time Markov chain.

Figure 5.1: Structure of an order fulfilment system in the exact modelling approach.

In contrast, the simplified modelling approach abstracts from the multi-stage character of order processing by combining the multiple processing steps of the order fulfilment system to one aggregated process (see Figure 5.2). Hence, each multi-stage order fulfilment system is modelled as a single-stage order fulfilment system consisting of one aggregated process, whose processing performance results from the processing performances of the different processing steps. Consequently, partially processed orders are neglected in the simplified model. An order is either completely unprocessed or completely processed. In this manner, this simplified modelling approach decreases both modelling complexity and computational effort of the discrete-time Markov chain.

Figure 5.2: Structure of an order fulfilment system in the simplified modelling approach.

## **5.4 Chapter Conclusion**

In this chapter, we selected and specified the modelling approach to model and analyse system behaviour and performance of a multi-stage, stochastic order fulfilment system with levelled order release and customer-required order deadlines. We chose a discrete-time Markov chain since discrete-time analytical models provide several advantages regarding accuracy, level of detail, and the description of real processes compared to continuous-time analytical models. Above all, a discrete-time Markov chain enables the exact calculation of complete probability distributions of stochastic performance measures. Thus, it provides a more detailed performance analysis of order fulfilment systems than continuous-time analytical models.

To model system behaviour of multi-stage, stochastic order fulfilment systems with levelled order release and customer-required order deadlines as a discretetime Markov chain, we consider the order types of the order fulfilment system separately. Moreover, we assume each order type to have a sequential processing sequence, and we limit the possible backlog duration of an order by the maximum accepted backlog duration. There are different modelling approaches to model the multi-stage character of order processing and the transmission of orders between successive processing steps in a multi-stage order fulfilment system as a discrete-time Markov chain. In this thesis, we differentiate between an exact and a simplified modelling approach. The exact and the simplified model for performance analysis will be introduced in Chapters 6 and 7. Furthermore, we will evaluate these models regarding modelling accuracy, accuracy of performance analysis, and memory and computation time requirements in Chapter 8.

# **6 Exact Model for Performance Analysis**

This chapter aims at introducing an exact model for performance analysis of multi-stage, stochastic order fulfilment systems with levelled order release and customer-required order deadlines. Section 6.1 provides several preliminary remarks about the exact model for performance analysis. We introduce the discrete-time Markov chain in Section 6.2. Based on the limiting distribution of the Markov chain, we derive exact formulas for various stochastic and deterministic, system- and customer-related performance measures of the order fulfilment system in Section 6.3. Section 6.4 discusses the memory and computation time requirements of the provided model. By summarising the results of this chapter, Section 6.5 provides an answer to the second research question of this thesis.

## **6.1 Preliminaries**

The exact model for performance analysis is based on the exact modelling approach. It considers a multi-stage, stochastic order fulfilment system with a finite set of processes P, a unique order type, and levelled order release considering a finite set of due dates K (see Figure 5.1). The order type and each process are specified by the parameters introduced in Sections 3.2.2 and 3.2.3, respectively.

We observe the order fulfilment system at discrete-time points in time t ∈ N<sup>0</sup> that are integer multiples of a constant time period. The time period corresponds to the scheduling interval of operational planning of the levelling concept (see Section 4.3). It commonly has a length of one day or one shift. In each time period, the order fulfilment system is affected by two stochastic processes:


### **6.1.1 Processing Performance per Time Period**

The processing performance per time period H<sup>p</sup> of process p ∈ P specifies the number of orders that can be completely processed at process p within one time period. At least (c<sup>p</sup> · lp,min) orders and at most (c<sup>p</sup> · lp,max) orders can be processed at process p within one time period. Hence, the range H<sup>p</sup> of the processing performance per time period H<sup>p</sup> is defined by

$$\mathcal{H}\_p = \left\{ (c\_p \cdot l\_{p,min}) , \dots , (c\_p \cdot l\_{p,max}) \right\}. \tag{6.1}$$

Assuming identical and independent distribution of the processing performance per time unit L<sup>p</sup> of process p, the probability distribution of the processing performance per time period H<sup>p</sup> of process p is computed as cp-fold convolution of the probability distribution of Lp.

We use the random vector

$$\mathbf{H} = \begin{pmatrix} H\_1 & \dots & H\_{p\_{max}} \end{pmatrix} \tag{6.2}$$

to describe the processing performance per time period of the order fulfilment system, whereby the random variable H<sup>p</sup> corresponds to the processing performance per time period of process p ∈ P. The range of H is given by

$$
\mathcal{H} = \mathcal{H}\_1 \times \dots \times \mathcal{H}\_{p\_{max}}.\tag{6.3}
$$

#### **6.1.2 Order Income per Time Period**

We specify the order income per time period using the random vector

$$\mathbf{G} = \begin{pmatrix} G\_{-R} & \dots & G\_{e\_{max}} \end{pmatrix},\tag{6.4}$$

whereby the random variable G<sup>k</sup> corresponds to the number of incoming orders per time period with a lead time of k ∈ K time periods. Its range G is defined based on the following conditions:


$$\mathcal{G} = \left\{ \mathbf{g} \in \mathbb{N}\_0^{(R + e\_{max} + 1)} \mid g\_k = 0 \right. \qquad \forall k \in \{-R, \dots, -1\} $$

$$\land g\_k \in \{0, \dots, a\_{max}\} \quad \forall k \in \{0, \dots, e\_{max}\} \qquad (6.5) $$

$$\land \sum\_{k \in \mathcal{K}} g\_k \in \mathcal{A} \right\}.$$

The probability P(G = g) depends on


$$P(\mathbf{G} = \mathbf{g}) = P\left(A = \sum\_{k \in \mathcal{K}} g\_k\right) \cdot \left[\prod\_{k=0}^{c\_{max}} \frac{(\sum\_{m=k}^{c\_{max}} g\_m)!}{g\_k! \cdot ((\sum\_{m=k}^{c\_{max}} g\_m) - g\_k)!}\right] \tag{6.6}$$

$$\cdot \left[\prod\_{k=0}^{c\_{max}} P(E = k)^{g\_k}\right] \qquad \forall \mathbf{g} \in \mathcal{G}.$$

#### **6.1.3 Stable Order Fulfilment System**

We assume the order fulfilment system to be stable. An order fulfilment system is stable if its order income-related utilisation U˜ is smaller than one. The order income-related utilisation U˜ is calculated as the ratio of the expected value of the number of incoming orders per time period E(A) to the smallest expected value of the processing performance per time period E(Hp) of the processes p ∈ P:

$$
\bar{U} = \frac{E(A)}{\min\_{p \in \mathcal{P}} E(H\_p)}.\tag{6.7}
$$

### **6.2 Discrete-time Markov Chain**

In the following, we introduce the discrete-time Markov chain that models system behaviour of a multi-stage, stochastic order fulfilment system with levelled order release and customer-required order deadlines according to the exact modelling approach.

#### **6.2.1 System State**

The system state X of the Markov chain depicts the number of unprocessed orders in the order fulfilment system at the beginning of any time period

$$\mathbf{X} = \begin{pmatrix} X\_{1, -R} & \dots & X\_{1, e\_{max}} \\ \vdots & \ddots & \vdots \\ X\_{p\_{max}, -R} & \dots & X\_{p\_{max}, e\_{max}} \end{pmatrix},\tag{6.8}$$

whereby Xp,k, p ∈ P, k ∈ K, specifies the number of unprocessed orders at process p with a due date of k time periods at the beginning of the time period.

#### **6.2.2 State Transition**

The state transition from an arbitrary state X<sup>t</sup> = x at the beginning of time period t to state Xt+1 = z at the beginning of time period (t + 1) consists of the following sub-steps:

	- Order processing at process p = 1,
	- . . . ,
	- Order processing at process pmax,

The sub-step of order processing is further subdivided into pmax sub-steps since we assume that order processing at the different processes is decoupled in time. Hence, order processing at process p does not start before order processing at the previous process (p − 1) is completed.

Based on these (pmax + 2) sub-steps, we explain the state transition in the following. We denote the interim states of the state transition by the variables y (m) , m ∈ {0, . . . ,(pmax + 2)}, whereby m is used as a counter. The initial value is given by y (0) = x. Note that y (p−1) specifies the number of unprocessed orders immediately before the start of order processing at process p ∈ P, whereas y (p) specifies the number of unprocessed orders after order processing at process p ∈ P is completed.

Sub-step of Order Processing The number of processed orders at process p in time period t depends on the realisation of the processing performance per time period H<sup>p</sup> = h<sup>p</sup> of process p in time period t. The number of processed orders with a due date of k time periods at process p in time period t

$$\min \left\{ y\_{p,k}^{(p-1)}; \max \left\{ 0; h\_p - \sum\_{l=-R}^{k-1} y\_{p,l}^{(p-1)} \right\} \right\}$$

equals either the number of unprocessed orders y (p−1) p,k immediately before the start of order processing at process p having a due date of k time periods or the residual processing performance remaining after all orders with a due date of (l < k) time periods have already been processed. Consequently, after order processing at process p is completed, the number of unprocessed orders at process p is given by

$$\begin{split} y\_{p,k}^{(p)} &= y\_{p,k}^{(p-1)} - \min \left\{ y\_{p,k}^{(p-1)}; \max \left\{ 0; h\_p - \sum\_{l=-R}^{k-1} y\_{p,l}^{(p-1)} \right\} \right\} \\ &= \max \left\{ 0; y\_{p,k}^{(p-1)} - \max \left\{ 0; h\_p - \sum\_{l=-R}^{k-1} y\_{p,l}^{(p-1)} \right\} \right\} \qquad \forall k \in \mathcal{K}, \end{split} \tag{6.9}$$

and the number of unprocessed orders at process (p + 1) equals

$$y\_{(p+1),k}^{(p)} = y\_{(p+1),k}^{(p-1)} + \min\left\{ y\_{p,k}^{(p-1)}; \max\left\{ 0; h\_p - \sum\_{l=-R}^{k-1} y\_{p,l}^{(p-1)} \right\} \right\} \tag{6.10}$$
 
$$\forall k \in \mathcal{K}.$$

Sub-step of Due Date Update Due to the transition from time period t to time period (t + 1), we have to update the due dates of all unprocessed orders at the end of time period t by reducing their due date by one time period:

$$\begin{aligned} y\_{p,k}^{(p\_{max}+1)} &= y\_{p,(k+1)}^{(p\_{max})} & \forall p \in \mathcal{P}, \forall k \in \mathcal{K} \\ y\_{p,e\_{max}}^{(p\_{max}+1)} &= 0 & \forall p \in \mathcal{P}. \end{aligned} \tag{6.11}$$

Consequently, unprocessed orders having a due date of (−R) time periods at the end of time period t become lost sales.

Sub-step of Order Income The incoming orders at the beginning of time period (t+ 1) are specified by the realisation of the order income per time period G = g at the beginning of time period (t+1). We add the incoming orders to the unprocessed orders at process p = 1 since the processes are ordered according to the processing sequence of the considered order type:

$$\begin{aligned} y\_{1,k}^{(p\_{max}+2)} &= y\_{1,k}^{(p\_{max}+1)} + g\_k & \forall k \in \mathcal{K} \\ y\_{p,k}^{(p\_{max}+2)} &= y\_{p,k}^{(p\_{max}+1)} & \forall p \in \mathcal{P} \end{aligned} \tag{6.12}$$

Transition Probability The transition probability results from the product of the probabilities of the aforementioned sub-steps since we assume that they are independent of each other. The sub-step of order processing depends on the realisation of the processing performance per time period H = h of the order fulfilment system in time period t. Thus, the probabilities P (p) , p ∈ P, are given as follows

$$P^{(p)} = P(H\_p = h\_p) \qquad \forall p \in \mathcal{P}.$$

The sub-step of due date update is deterministic,

$$P^{(p\_{max}+1)} = 1,$$

whereas the sub-step of order income depends on the realisation of the order income per time period G = g at the beginning of time period (t + 1):

$$P^{(p\_{max}+2)} = P(\mathbf{G} = \mathbf{g}).$$

By combining these probabilities, we define the transition probability of the state transition from an arbitrary state X<sup>t</sup> = x at the beginning of time period t to state X<sup>t</sup>+1 = z at the beginning of time period (t + 1) as follows

$$P(\mathbf{X}^{t+1} = \mathbf{z} \mid \mathbf{X}^t = \mathbf{x}) = \sum\_{\{\mathbf{g}, \mathbf{h}\} \in \mathcal{T}(\mathbf{x}, \mathbf{z})} \left[ \prod\_{p \in \mathcal{P}} P(H\_p = h\_p) \right] \cdot P(\mathbf{G} = \mathbf{g}) \tag{6.13}$$
 
$$\forall \mathbf{x}, \mathbf{z} \in \mathcal{X},$$

whereby set I (x, z) contains all tuples (g, h) ∈ G × H resulting in a state transition from X<sup>t</sup> = x to X<sup>t</sup>+1 = z.

#### **6.2.3 State Space**

The state space X of the Markov chain is non-negative as it is impossible to observe a negative number of orders. The upper bound of the state space

$$\mathbf{O} = \begin{pmatrix} O\_{1,-R} & \dots & O\_{1,e\_{max}} \\ \vdots & \ddots & \vdots \\ O\_{p\_{max},-R} & \dots & O\_{p\_{max},e\_{max}} \end{pmatrix} \\ \tag{6.14}$$

results from the structure of the state transition. We observe the maximum possible number of unprocessed orders Op,k at process p ∈ P with a due date of k ∈ K time periods in any time period if the following conditions hold:


The maximum possible number of incoming orders at process p = 1 equals the maximum possible number of incoming orders per time period amax for non-negative due dates and zero for negative due dates. At all subsequent processes p ∈ P \ {1}, it corresponds to the maximum possible processing performance per time period h(p−1),max of the previous processing step (p−1). The minimum possible residual processing performance at process p equals the minimum possible processing performance per time period hp,min of process p for a due date of (−R) time periods. For due dates of (k > −R) time periods, it is zero. The maximum possible number of orders remaining from previous time periods at process p is zero if the due date corresponds to maximum possible lead time of emax time periods. For due dates of (k < emax) time periods, it corresponds to the maximum possible number of unprocessed orders with a due date of (k + 1) time periods.

Thus, the components Op,k of the upper bound O are calculated as follows

O1,k = (emax + 1) · amax ∀k ∈ {−R, . . . , −1} O1,k = (emax − k + 1) · amax ∀k ∈ {0, . . . , emax} Op,k = (emax + |k|) · h(p−1),max ∀p ∈ P \ {1}, ∀k ∈ {−R, . . . , −1} Op,k = (emax − k + 1) · h(p−1),max ∀p ∈ P \ {1}, ∀k ∈ {0, . . . , emax} . (6.15)

#### **6.2.4 Limiting Distribution**

The Markov chain has a finite state space (see Section 6.2.3), and it is possible to reach every state of the Markov chain from every other state either by a singlestep transition or by an indirect transition via a finite number of other states. In the case of unreachable states, these states are excluded from the computations (see Section 10.1.1). We then find an irreducible subset of the state space. We assume that there is a tuple (g, h) ∈ G × H that meets the following conditions:

• The processing performance per time period is the same at every process p ∈ P:

$$h\_p = h \qquad \forall p \in \mathcal{P}.$$

• The total number of incoming orders per time period equals the processing performance per time period:

$$\sum\_{k \in \mathcal{K}} g\_k = h.$$

Thus, the irreducible subset of the state space contains state x with

$$\begin{aligned} x\_{1,k} &= g\_k \quad \forall k \in \mathcal{K} \\ x\_{p,k} &= 0 \quad \forall p \in \mathcal{P} \\ \sum\_{k \in \mathcal{K}} x\_{1,k} &= h, \end{aligned}$$

for which

$$P(\mathbf{X}^{t+1} = \mathbf{x} \mid \mathbf{X}^t = \mathbf{x}) = \left[ \prod\_{p \in \mathcal{P}} P(H\_p = h) \right] \cdot P(\mathbf{G} = \mathbf{g}) > 0$$

holds. Consequently, state x as well as all other states of the irreducible subset of the state space are aperiodic. In conclusion, the limiting distribution π˜ of the Markov chain exists, and we calculate the limiting distribution π˜ of a given initial state as described in Section 5.2.2.

## **6.3 Performance Measures of the Order Fulfilment System**

In the following, we derive several stochastic and deterministic, system- and customer-related performance measures of the order fulfilment system from the limiting distribution of the Markov chain. Table 6.1 provides an overview of the calculated performance measures of the order fulfilment system. The probability of being in state X = x after a sufficiently long time, in steady-state, is denoted by P(X = x) := π˜x.

#### **6.3.1 Number of Unprocessed Orders**

For each process p ∈ P, we differentiate between


Since order processing at the different processing steps of the order fulfilment system is decoupled in time, the values of the performances measures Q<sup>p</sup> and Q˜ p only equal for the first process. Q<sup>p</sup> is a suitable performance measure to quantify the overall order backlog in the order fulfilment system at any point in time, and Q˜ <sup>p</sup> is a meaningful performance measure for buffer sizing at the processing steps p ∈ P of the system.


Table 6.1: Performance measures of the order fulfilment system.

The number of unprocessed orders Q<sup>p</sup> at process p at the beginning of a time period is the sum of the number of unprocessed orders xp,k at process p with a due date of k time periods for all due dates k ∈ K. Its range Q<sup>p</sup> is limited downwards by zero and upwards by the sum of the upper limit of the state space Op,k at process p with a due date of k time periods for all due dates k ∈ K:

$$\mathcal{Q}\_p = \left\{ 0, \ldots, \left( \sum\_{k \in \mathcal{K}} O\_{p,k} \right) \right\}.\tag{6.16}$$

The probability distribution of Qp, p ∈ P, is calculated as follows

$$\begin{aligned} P(Q\_p = q\_p) &= \sum\_{\mathbf{x} \in \mathcal{Z}\_{\mathcal{Q}}(p, q\_p)} P(\mathbf{X} = \mathbf{x}) & \forall q\_p \in \mathcal{Q}\_p\\ \mathcal{Z}\_{\mathcal{Q}}(p, q\_p) &= \left\{ \mathbf{x} \in \mathcal{X} \mid \sum\_{k \in \mathcal{K}} x\_{p, k} = q\_p \right\}. \end{aligned} \tag{6.17}$$

Furthermore, the number of unprocessed orders Q in the order fulfilment system at the beginning of a time period is the sum of the number of unprocessed orders xp,k at process p with a due date of k time periods for all processes p ∈ P and all due dates k ∈ K. Its range Q is limited downwards by zero and upwards by the sum of the upper limit of the state space Op,k at process p with a due date of k time periods for all processes p ∈ P and all due dates k ∈ K:

$$\mathcal{Q} = \left\{ 0, \ldots, \left. \left( \sum\_{p \in \mathcal{P}} \sum\_{k \in \mathcal{K}} O\_{p,k} \right) \right\} \right. \tag{6.18}$$

The probability distribution of Q is calculated as follows

$$\begin{aligned} P(Q = q) &= \sum\_{\mathbf{x} \in \mathcal{X}\_{\mathcal{Q}}(q)} P(\mathbf{X} = \mathbf{x}) & \forall q \in \mathcal{Q} \\ \mathcal{Z}\_{\mathcal{Q}}(q) &= \left\{ \mathbf{x} \in \mathcal{X} \mid \sum\_{p \in \mathcal{P}} \sum\_{k \in \mathcal{K}} x\_{p,k} = q \right\}. \end{aligned} \tag{6.19}$$

The number of unprocessed orders Q˜ <sup>p</sup> at process p immediately before the start of order processing at process p corresponds to the sum of unprocessed orders y (p−1) p,k at process p with a due date of k time periods immediately before the start of order processing at process p for all due dates k ∈ K. Its range Q˜ <sup>p</sup> is limited downwards by zero. At process p = 1, the upper bound of Q˜ <sup>1</sup> is given by the sum of the upper limit of the state state O1,k at process p = 1 with a due date of k time periods for all due dates k ∈ K. Otherwise, the upper bound of Q˜ <sup>p</sup> additionally incorporates the maximum possible processing performance per time period h(p−1),max of the previous processing step (p − 1):

$$\bar{\mathcal{Q}}\_p = \begin{cases} \{0, \dots, \left(\sum\_{k \in \mathcal{K}} O\_{p,k}\right)\} & p = 1\\ \{0, \dots, \left(\sum\_{k \in \mathcal{K}} O\_{p,k} + h\_{(p-1),max}\right)\} & p \in \mathcal{P} \nmid \{1\}. \end{cases} \tag{6.20}$$

The probability distribution of Q˜ <sup>p</sup>, p ∈ P, is calculated as follows

$$\begin{split} P(\bar{Q}\_{p} = \bar{q}\_{p}) &= \sum\_{\{\mathbf{x}, \mathbf{h}\} \in \mathcal{Z}\_{\mathcal{Q}}(p, \bar{q}\_{p})} P(\mathbf{X} = \mathbf{x}) \cdot P(\mathbf{H} = \mathbf{h}) \qquad \forall \bar{q}\_{p} \in \bar{\mathcal{Q}}\_{p} \\ \mathcal{Z}\_{\mathcal{Q}}(p, \bar{q}\_{p}) &= \left\{ (\mathbf{x}, \mathbf{h}) \in \mathcal{X} \times \mathcal{H} \mid \sum\_{k \in \mathcal{K}} y\_{p,k}^{(p-1)} = \bar{q}\_{p} \right\}. \end{split} \tag{6.21}$$

### **6.3.2 Number of Lost Sales**

Lost sales occur in the order fulfilment system since the set of due dates K is limited downwards by the maximum backlog duration R. Unprocessed orders at process p with a due date of (−R) time periods become lost sales in any time period t if their number exceeds the processing performance per time period of process p in time period t. Thus, some of these orders are still unprocessed at the end of time period t and their due dates fall below the minimum possible due date when the due dates are updated at the end of time period t. These unprocessed orders are removed from the order backlog of the order fulfilment system without being processed. The performance measures Sp, p ∈ P, and S quantify the number of lost sales per time period at process p and the total number of lost sales per time period in the order fulfilment system, respectively.

The number of lost sales per time period S<sup>p</sup> at process p corresponds to the number of unprocessed orders y (p) p,−<sup>R</sup> with a due date of (−R) time periods remaining after order processing at process p is completed. Its range S<sup>p</sup> is limited by zero and the upper limit of the state space Op,−<sup>R</sup> at process p with a due date of (−R) time periods:

$$\mathcal{S}\_p = \{0, \dots, O\_{p,-R}\}\,. \tag{6.22}$$

The probability distribution of Sp, p ∈ P, is calculated as follows

$$\begin{split} P(S\_p = s\_p) &= \sum\_{\mathbf{(x,h)} \in \mathcal{T}\_S(p, s\_p)} P(\mathbf{X} = \mathbf{x}) \cdot P(\mathbf{H} = \mathbf{h}) \quad \forall s\_p \in \mathcal{S}\_p \\ \mathcal{Z}\_{\mathbf{S}}(p, s\_p) &= \left\{ (\mathbf{x}, \mathbf{h}) \in \mathcal{X} \times \mathcal{H} \mid y\_{p, -R}^{\{p\}} = s\_p \right\}. \end{split} \tag{6.23}$$

The total number of lost sales per time period S in the order fulfilment system is the sum of the number of unprocessed orders y (p) p,−<sup>R</sup> with a due date of (−R) time periods remaining after order processing is completed for all processes p ∈ P. Its range S is limited downwards by zero and upwards by the sum of the upper limit of the state space Op,−<sup>R</sup> at process p with a due date of (−R) time periods for all processes p ∈ P:

$$\mathcal{S} = \left\{ 0, \ldots, \sum\_{p \in \mathcal{P}} O\_{p, -R} \right\}. \tag{6.24}$$

The probability distribution of S is calculated as follows

$$\begin{split}P(S=s) &= \sum\_{\{\mathbf{x}, \mathbf{h}\} \in \mathcal{T}\_{\mathcal{S}}(s)} P(\mathbf{X}=\mathbf{x}) \cdot P(\mathbf{H}=\mathbf{h}) \qquad \forall s \in \mathcal{S} \\ &\mathcal{T}\_{\mathcal{S}}(s) = \left\{ (\mathbf{x}, \mathbf{h}) \in \mathcal{X} \times \mathcal{H} \, \big|\, \sum\_{p \in \mathcal{P}} y\_{p,-R}^{(p)} = s \right\}. \end{split} \tag{6.25}$$

#### **6.3.3 Performance Balance**

The performance balance W<sup>p</sup> of process p ∈ P specifies the difference between the processing performance per time period of process p and the number of unprocessed orders immediately before the start of order processing at process p. It is measured in number of unprocessed orders, and it is a suitable performance measure to evaluate the workload at process p. A positive performance balance at process p ∈ P indicates that there is some idle time at process p ∈ P since the processing performance provided per time period of process p exceeds the number of unprocessed orders per time period immediately before the start of order processing at process p. Otherwise, in the case of a negative performance balance, the provided processing performance at process p ∈ P is fully utilised since the number of unprocessed orders per time period immediately before the start of order processing at process p exceeds the processing performance per time period of process p.

We observe a performance balance of w<sup>p</sup> orders at process p if the processing performance per time period h<sup>p</sup> of process p differs by w<sup>p</sup> orders from the sum of the number of unprocessed orders y (p−1) p,k at process p with a due date of k time periods immediately before the start of order processing at process p for all due dates k ∈ K. The range W<sup>p</sup> of the performance balance of process p is limited downwards by the negated maximum possible number of unprocessed orders at process p immediately before the start of order processing at process p. Its upper bound is given by the maximum possible processing performance per time period of process p. Hence, the range W<sup>p</sup> is defined as follows

$$\mathcal{W}\_p = \begin{cases} \left\{ -\left( \sum\_{k \in \mathcal{K}} O\_{p,k} \right), \dots, h\_{p,max} \right\} & p = 1\\ \left\{ -\left( \sum\_{k \in \mathcal{K}} O\_{p,k} + h\_{\{p-1\},max} \right), \dots, h\_{p,max} \right\} & p \in \mathcal{P} \nmid \{1\}. \end{cases} \tag{6.26}$$

The probability distribution of Wp, p ∈ P, is computed as follows

$$\begin{split} P(W\_{p} = w\_{p}) &= \sum\_{\{\mathbf{x}, \mathbf{h}\} \in \mathcal{Z}\_{\mathcal{W}}(p, w\_{p})} P(\mathbf{X} = \mathbf{x}) \cdot P(\mathbf{H} = \mathbf{h}) \qquad \forall w\_{p} \in \mathcal{W}\_{p} \\ \mathcal{Z}\_{\mathcal{W}}(p, w\_{p}) &= \left\{ (\mathbf{x}, \mathbf{h}) \in \mathcal{X} \times \mathcal{H} \mid h\_{p} - \sum\_{k \in \mathcal{K}} y\_{p,k}^{\{p-1\}} = w\_{p} \right\}. \end{split} \tag{6.27}$$

#### **6.3.4 Utilisation**

The order backlog-related utilisation U<sup>p</sup> of process p describes the relation between the processing performance per time period and the number of unprocessed orders immediately before the start of order processing at process p. It considers the ratio of the number of unprocessed orders immediately before the start of order processing at process p to the processing performance per time period of process p. Its range is limited to the interval [0,1]. We calculate the order backlog-related utilisation of process p as follows

$$U\_p = \sum\_{(\mathbf{x}, \mathbf{h}) \in \mathcal{X} \times \mathcal{H}} \min \left\{ 1; \frac{\sum\_{k \in \mathcal{K}} y\_{p,k}^{(p-1)}}{h\_p} \right\} \cdot P(\mathbf{X} = \mathbf{x}) \cdot P(\mathbf{H} = \mathbf{h}).\tag{6.28}$$

#### **6.3.5 Number of Processed Orders**

The number of processed orders per time period is a suitable performance measure to quantify the throughput of any process p ∈ P (performance measure Fp) as well as the throughput of the whole system (performance measure F).

We initially calculate the processed orders per time period

$$\mathbf{M} = \begin{pmatrix} M\_{1, -R} & \dots & M\_{1, e\_{max}} \\ \vdots & \ddots & \vdots \\ M\_{p\_{max}, -R} & \dots & M\_{p\_{max}, e\_{max}} \end{pmatrix},\tag{6.29}$$

whereby Mp,k specifies the number of processed orders per time period at process p with a due date of k time periods at their time of processing. Its range M is defined based on the following conditions:


M = ( m ∈ N pmax <sup>0</sup> ×N (R+emax+1) 0 | mp,k ∈ {0, . . . , hp,max} ∀p ∈ P, ∀k ∈ K ∧ X k∈K <sup>m</sup>p,k ∈ H<sup>p</sup> <sup>∀</sup><sup>p</sup> ∈ P) . (6.30)

The number of processed orders per time period at process p with a due date of k time periods at their time of processing

$$\min \left\{ y\_{p,k}^{(p-1)}; \max \left\{ 0; h\_p - \sum\_{l=-R}^{k-1} y\_{p,l}^{(p-1)} \right\} \right\}$$

equals either the number of unprocessed orders y (p−1) p,k at process p with a due date of k time periods immediately before the start of order processing at process p or the residual processing performance of process p remaining after all orders with a due date of (l < k) time periods have already been processed.

The probability distribution of M is calculated as follows

P(M = m) = X (x,h)∈IM(m) P(X = x) · P(H = h) ∀m ∈ M <sup>I</sup>M(m) = ( (x, h) ∈ X × H | min ( y (p−1) p,k ; max ( 0; h<sup>p</sup> − kX−1 l=−R y (p−1) p,l )) <sup>=</sup> <sup>m</sup>p,k <sup>∀</sup><sup>p</sup> ∈ P, <sup>∀</sup><sup>k</sup> ∈ K) . (6.31)

The number of processed orders per time period F<sup>p</sup> at process p corresponds to the sum of the number of processed orders mp,k at process p with a due date of k time periods at their time of processing for all due dates k ∈ K. Its range F<sup>p</sup> is limited downwards by zero and upwards by the maximum possible processing performance per time period hp,max of process p:

$$\mathcal{F}\_p = \{0, \dots, h\_{p,max}\}.\tag{6.32}$$

The probability distribution of Fp, p ∈ P, is computed as follows

$$\begin{aligned} P(F\_p = f\_p) &= \sum\_{\mathbf{m} \in \mathbb{Z}\_F(p, f\_p)} P(\mathbf{M} = \mathbf{m}) & \forall f\_p \in \mathcal{F}\_p \\ \mathcal{Z}\_{\mathcal{F}}(p, f\_p) &= \left\{ \mathbf{m} \in \mathcal{M} \mid \sum\_{k \in \mathcal{K}} m\_{p,k} = f\_p \right\}. \end{aligned} \tag{6.33}$$

The number of completed orders per time period F corresponds to the number of processed orders per time period at process pmax of the order fulfilment system. Hence, range and probability distribution of F correspond to the ones of Fpmax :

$$\mathcal{F} = \mathcal{F}\_{p\_{max}} \tag{6.34}$$

$$P(F = f) = P\left(F\_{p\_{max}} = f\right) \qquad \forall f \in \mathcal{F}.\tag{6.35}$$

### **6.3.6 Time Difference to Order Deadline**

The time difference to order deadline D specifies the time difference between the deadline and the actual time of completion of an order. A negative time difference to order deadline indicates that the order is completed after its deadline, whereas a positive time difference to order deadline indicates that the order is completed before reaching its deadline. The range D of the time difference to order deadline corresponds to the set of due dates:

$$\mathcal{D} = \mathcal{K}.\tag{6.36}$$

The probability P(D = d) that a completed order has a time difference of d time periods to its deadline at its time of completion is proportional to the weighted sum of the number of processed orders per time period mpmax,d at process pmax with a due date of d time periods at their time of completion for all realisations m ∈ M. By normalising this weighted sum, we calculate the probability distribution of D as follows

$$P(D=d) = \frac{\sum\_{\mathbf{m}\in\mathcal{M}} m\_{p\_{max},d} \cdot P(\mathbf{M}=\mathbf{m})}{\sum\_{\mathbf{m}\in\mathcal{M}} \left(\sum\_{k\in\mathcal{K}} m\_{p\_{max},k}\right) \cdot P(\mathbf{M}=\mathbf{m})} \qquad \forall d \in \mathcal{D}.\tag{6.37}$$

Based on the time difference to order deadline D, we specify the performance measures backlog duration Dbacklog and time buffer Dbuf f er of a completed order in the following.

The backlog duration of a completed order Dbacklog equals the number of time periods by which its time of completion exceeds its deadline. Its range Dbacklog is limited downwards by one and upwards by the maximum backlog duration R:

$$\mathcal{D}^{backlog} = \{1, \ldots, R\}.\tag{6.38}$$

The probability P(Dbacklog = d) that a completed order has a backlog duration of d time periods at its time of completion is proportional to the probability P(D = −d) that a completed order has a time difference to deadline of (−d) time periods. By normalising the probabilities P(D = −d), d ∈ Dbacklog, the probability distribution of Dbacklog is computed as follows

$$P(D^{backlog} = d) = \frac{P(D = -d)}{\sum\_{k=-R}^{-1} P(D = k)} \qquad \forall d \in \mathcal{D}^{backlog}.\tag{6.39}$$

The time buffer of a completed order Dbuf f er corresponds to the number of time periods that remain at its time of completion until reaching its deadline. Its range Dbuf f er corresponds to the subset of non-negative due dates:

$$\mathcal{D}^{buffer} = \{0, \ldots, e\_{max}\}\,. \tag{6.40}$$

The probability P(Dbuf f er = d) that a completed order has a time buffer of d time periods at its time of completion is proportional to the probability P(D = d) that a completed order has a time difference to deadline of d time periods. By normalising the probabilities P(D = d), d ∈ Dbuf f er, the probability distribution of Dbuf f er is computed as follows

$$P(D^{buffer} = d) = \frac{P(D = d)}{\sum\_{k=0}^{c\_{max}} P(D = k)} \qquad \forall d \in \mathcal{D}^{buffer}.\tag{6.41}$$

#### **6.3.7 Service Level**

The service level of an order fulfilment system specifies the proportion of ontime completed orders on the total number of orders. Its range is limited to the interval [0,1]. Regarding the level of detail, we differentiate three types of service level:


• γ-service level SL<sup>γ</sup> considers the proportion of backorders as well as their backlog duration.

We define α-service level SL<sup>α</sup> as the probability that none of the completed orders per time period has a backlog duration. Thus, we consider the probability that the sum of processed orders per time period mpmax,k at process pmax with a negative due date k ∈ {−R, . . . , −1} at its time of completion equals zero:

$$\begin{aligned} SL\_{\alpha} &= \sum\_{\mathbf{m} \in \mathbb{Z}\_{\alpha}} P(\mathbf{M} = \mathbf{m}) \\ \mathcal{I}\_{\alpha} &= \left\{ \mathbf{m} \in \mathcal{M} \mid \sum\_{k=-R}^{-1} m\_{p\_{max},k} = 0 \right\}. \end{aligned} \tag{6.42}$$

β-service level SL<sup>β</sup> is the proportion of on-time completed orders on the total number of outgoing orders. The number of on-time completed orders per time period corresponds to the number of processed orders per time period at process pmax with a non-negative due date k ∈ {0, . . . , emax}. The total number of outgoing orders per time period is the sum of the total number of processed orders per time period at process pmax and the total number of lost sales per time period. We calculate SL<sup>β</sup> as follows

$$\begin{split} SL\_{\beta} &= \sum\_{\mathbf{x} \in \mathcal{X}} \sum\_{\mathbf{h} \in \mathcal{H}} \frac{\Phi\_{\beta}(\mathbf{x}, \mathbf{h})}{\Psi\_{\beta}(\mathbf{x}, \mathbf{h})} \cdot P(\mathbf{X} = \mathbf{x}) \cdot P(\mathbf{H} = \mathbf{h}) \\ \Phi\_{\beta}(\mathbf{x}, \mathbf{h}) &= \sum\_{k=0}^{c\_{\max}} \min \left\{ y\_{p\_{\max},k}^{(p\_{\max}-1)}; \max \left\{ 0; h\_{p\_{\max}} - \sum\_{l=-R}^{k-1} y\_{p\_{\max},l}^{(p\_{\max}-1)} \right\} \right\} \\ \Psi\_{\beta}(\mathbf{x}, \mathbf{h}) &= \sum\_{k \in \mathcal{K}} \min \left\{ y\_{p\_{\max},k}^{(p\_{\max}-1)}; \max \left\{ 0; h\_{p\_{\max}} - \sum\_{l=-R}^{k-1} y\_{p\_{\max},l}^{(p\_{\max}-1)} \right\} \right\} \\ &\quad + \sum\_{p \in \mathcal{P}} y\_{p,-R}^{(p)} .\end{split} \tag{6.43}$$

We define γ-service level SL<sup>γ</sup> as complementary of the proportion of completed orders having a backlog duration and the total number of lost sales on the total number of outgoing orders, whereby each component of the formula is weighted by its corresponding backlog duration. The number of completed orders per time period having a backlog duration corresponds to the number of processed orders per time period at process pmax with a negative due date k ∈ {−R, . . . , −1}. It is weighted by its backlog duration of |k| time periods. Furthermore, we define the weight of the total number of lost sales per time period by a backlog duration of (R + 1) time periods. The total number of outgoing orders per time period consists of the total number of processed orders at process pmax, which is weighted by the maximum backlog duration of R time periods, and the total number of lost sales per time period, weighted by a backlog duration of (R + 1) time periods. We calculate SL<sup>γ</sup> as follows

$$\begin{split} SL\_{\gamma} &= 1 - \left( \sum\_{\mathbf{x} \in \mathcal{X}} \sum\_{\mathbf{h} \in \mathcal{H}} \frac{\Phi\_{\gamma}(\mathbf{x}, \mathbf{h})}{\Psi\_{\gamma}(\mathbf{x}, \mathbf{h})} \cdot P(\mathbf{X} = \mathbf{x}) \cdot P(\mathbf{H} = \mathbf{h}) \right) \\ \Phi\_{\gamma}(\mathbf{x}, \mathbf{h}) &= \sum\_{k=-R}^{-1} |k| \cdot \min \left\{ y\_{p\_{max},k}^{(p\_{max}-1)}; \max \left\{ 0; h\_{p\_{max}} - \sum\_{l=-R}^{k-1} y\_{p\_{max},l}^{(p\_{max}-1)} \right\} \right\} \\ &\quad + (R+1) \cdot \sum\_{p \in \mathcal{P}} y\_{p,-R}^{(p)} \end{split}$$
 
$$\begin{split} \Psi\_{\gamma}(\mathbf{x}, \mathbf{h}) &= R \cdot \sum\_{k \in \mathcal{K}} \min \left\{ y\_{p\_{max},k}^{(p\_{max}-1)}; \max \left\{ 0; h\_{p\_{max}} - \sum\_{l=-R}^{k-1} y\_{p\_{max},l}^{(p\_{max}-1)} \right\} \right\} \\ &\quad + (R+1) \cdot \sum\_{p \in \mathcal{P}} y\_{p,-R}^{(p)} . \end{split} \tag{6.44}$$

To verify the implementation of the exact model for performance analysis, we conduct a model comparison between the exact model and a simulation model. The simulation model depicts system behaviour of a multi-stage, stochastic order fulfilment system with levelled order release and customer-required deadlines, analogous to the Markov chain. Based on the simulation results, we calculate the same performance measures of the order fulfilment system as in the exact model (see Table 6.1). A detailed description of the simulation model is provided in Appendix B. The results of this model comparison confirm that the exact model is implemented correctly. Details on the design and the results of the model comparison are given in Appendix D.

## **6.4 Memory and Computation Time Requirements**

Computing the exact model for performance analysis consists of calculating (1) the state space, (2) the transition matrix, (3) the limiting distribution of the Markov chain, and (4) the performance measures of the order fulfilment system. The memory and computation time requirements of these steps are predominantly driven by the size of the state space of the Markov chain. As the state space increases, the transition matrix becomes larger. Thus, larger sets of linear equations have to be solved to obtain the limiting distribution. Finally, calculating the performance measures of the order fulfilment system requires more memory and computation time since they are derived from the limiting distribution.

The size of the state space |X | of the Markov chain is derived from the lower and upper limit of the state space X (see Section 6.2.3) as follows

$$\begin{split} |\mathcal{X}| &= \prod\_{p \in \mathcal{P}} \prod\_{k \in \mathcal{K}} \left( O\_{p,k} + 1 \right) \\ |\mathcal{X}| &= \left[ \left( e\_{max} + 1 \right) \cdot a\_{max} + 1 \right]^R \cdot \prod\_{k=0}^{e\_{max}} \left( \left( e\_{max} - k + 1 \right) \cdot a\_{max} + 1 \right) \\ &\cdot \prod\_{p=2}^{p\_{max}} \left[ \prod\_{k=-R}^{-1} \left( \left( e\_{max} + |k| \right) \cdot h\_{(p-1),max} + 1 \right) \\ &\cdot \prod\_{k=0}^{e\_{max}} \left( \left( e\_{max} - k + 1 \right) \cdot h\_{(p-1),max} + 1 \right) \right]. \end{split} \tag{6.45}$$

An upper limit of the size of state space is obtained by the following estimations

$$\begin{aligned} |\mathcal{X}| &\le \left[ (e\_{\max} + 1) \cdot a\_{\max} + 1 \right]^{\{R + e\_{\max} + 1\}} \\ &\cdot \prod\_{p=2}^{p\_{\max}} \left[ \left( (e\_{\max} + R) \cdot h\_{\{p-1\}, \max} + 1 \right)^R \\ &\cdot \left( (e\_{\max} + 1) \cdot h\_{\{p-1\}, \max} + 1 \right)^{\{e\_{\max} + 1\}} \right] \end{aligned} \tag{6.46}$$

$$|\mathcal{X}| \le \left[ (e\_{max} + 1) \cdot a\_{max} + 1 \right]^{(R + e\_{max} + 1)}$$

$$\cdot \left[ (e\_{max} + R) \cdot h\_{max} + 1 \right]^{R \cdot (p\_{max} - 1)} \tag{6.47}$$

$$\cdot \left[ (e\_{max} + 1) \cdot h\_{max} + 1 \right]^{(e\_{max} + 1)(p\_{max} - 1)}$$

with

$$h\_{max} = \max\_{p=1,\dots,(p\_{max}-1)} \{h\_{p,max}\}\dots$$

From equation (6.47), we recognize that the maximum backlog duration R, the number of processes pmax, and the maximum possible values of the number of incoming orders per time period A, the lead time E, and the processing performances per time period H<sup>p</sup> of the processes p ∈ P determine the size of the state space of the Markov chain. The maximum backlog duration R, the maximum possible lead time emax, and the number of processes pmax have a major impact on the size of the state space since they determine the exponents in equation (6.47).

Apart from the size of the state space, the computational effort of calculating the transition matrix depends on the discretisation of the stochastic parameters. The discretisation of a stochastic parameter quantifies the number of its possible realisations. It results from the lower and upper limit of the range of the stochastic parameter. The higher the discretisation of the stochastic parameters G and H of the Markov chain, which results from the discretisation of the stochastic parameters A, E, and Lp, p ∈ P, of the order fulfilment system, the more tuples (g, h) of possible realisations exist. Each of these tuples has to be considered when calculating the entries of the transition matrix (see equation (6.13)).

## **6.5 Chapter Conclusion**

In this chapter, we introduced an exact analytical model for performance analysis of multi-stage, stochastic order fulfilment systems with levelled order release and customer-required order deadlines by


Based on the calculated performance measures (see Table 6.1), such as system throughput, service level, utilisation, and backlog duration and time buffer of a completed order, the developed model enables an exact and comprehensive performance analysis of multi-stage, stochastic order fulfilment systems with levelled order release and customer-required order deadlines. Consequently, the exact model for performance analysis provides an answer to the second research question of the thesis:

How can we determine the performance of multi-stage, stochastic order fulfilment systems with levelled order release and customerrequired order deadlines?

# **7 Simplified Model for Performance Analysis**

This chapter aims at introducing a simplified model for performance analysis of multi-stage, stochastic order fulfilment systems with levelled order release and customer-required order deadlines. Section 7.1 provides several preliminary remarks about the simplified model for performance analysis. We introduce the discrete-time Markov chain in Section 7.2 and calculate multiple performance measures of the order fulfilment system in Section 7.3. Section 7.4 discusses the memory and computation time requirements of the provided model and Section 7.5 summarises the results of this chapter.

## **7.1 Preliminaries**

The simplified model for performance analysis is based on the simplified modelling approach. It considers a multi-stage, stochastic order fulfilment system with a finite set of processes P, a unique order type, and levelled order release considering a finite set of due dates K. The processing steps p ∈ P are combined to one aggregated process (see Figure 5.2). The order type and each process is specified by the parameters introduced in Sections 3.2.2 and 3.2.3.

We observe the order fulfilment system at discrete-time points in time t ∈ N<sup>0</sup> that are integer multiples of a constant time period. The time period corresponds to the scheduling interval of operational planning of the levelling concept (see Section 4.3). It commonly has a length of one day or one shift. In each time period, the order fulfilment system is affected by two stochastic processes:


### **7.1.1 Aggregated Processing Performance per Time Period**

The aggregated processing performance per time period of the order fulfilment system specifies the number of orders that can be completely processed in the order fulfilment system within one time period. We calculate the aggregated processing performance per time period H˜ as the minimum of the processing performances per time period of the processes p ∈ P:

$$
\bar{H} = \min\_{p \in \mathcal{P}} \left\{ H\_p \right\}. \tag{7.1}
$$

Hence, its range H˜ results from the smallest value of the minimum possible processing performances per time period hp,min of the processes p ∈ P and the smallest value of the maximum possible processing performances per time period hp,max of the processes p ∈ P as follows

$$\bar{\mathcal{H}} = \left\{ \min\_{p \in \mathcal{P}} \{ h\_{p,min} \} , \dots , \min\_{p \in \mathcal{P}} \{ h\_{p,max} \} \right\}. \tag{7.2}$$

Assuming that the processing performance per time period of the processes p ∈ P are independent of each other, the probability distribution of H˜ is given as follows

$$\begin{split}P(\bar{H}=\bar{h}) &= \sum\_{\mathbf{h}\in\mathcal{Z}\_{\tilde{\mathcal{H}}}(\bar{h})} \left[\prod\_{p\in\mathcal{P}} P(H\_p = h\_p)\right] \quad \forall \bar{h}\in\bar{\mathcal{H}}\\ \mathcal{Z}\_{\tilde{\mathcal{H}}}(\bar{h}) &= \left\{\mathbf{h}\in\mathcal{H} \,|\,\min\_{p\in\mathcal{P}}\{h\_p\}=\bar{h}\right\}.\end{split} \tag{7.3}$$

#### **7.1.2 Order Income per Time Period**

The order income per time period is specified by the random vector G, analogous to the exact model (see Section 6.1.2). The range and the probability distribution of G are given by equations (6.5) and (6.6), respectively.

### **7.1.3 Stable Order Fulfilment System**

We assume the order fulfilment system to be stable. An order fulfilment system is stable, if its order income-related utilisation U˜ is smaller than one. In the simplified model, the order income-related utilisation U˜ is calculated as the ratio of the expected value of the number of incoming orders per time period E(A) to the expected value of the aggregated processing performance per time period E(H˜ ):

$$
\bar{U} = \frac{E(A)}{E(\bar{H})}.\tag{7.4}
$$

### **7.2 Discrete-time Markov Chain**

In the following, we introduce the discrete-time Markov chain modelling system behaviour of a multi-stage, stochastic order fulfilment system with levelled order release and customer-required deadlines according to the simplified modelling approach.

#### **7.2.1 System State**

The system state X of the Markov chain depicts the number of unprocessed orders in the order fulfilment system at the beginning of any time period. Since order fulfilment systems are modelled as single-stage systems in the simplified model, we neglect the process index p. Hence, the definition of system state simplifies to

$$\mathbf{X} = \begin{pmatrix} X\_{-R} & \dots & X\_{e\_{max}} \end{pmatrix},\tag{7.5}$$

whereby Xk, k ∈ K, specifies the number of unprocessed orders with a due date of k time periods at the beginning of the time period.

#### **7.2.2 State Transition**

The state transition from an arbitrary state X<sup>t</sup> = x at the beginning of time period t to state Xt+1 = z at the beginning of time period (t + 1) consists of the same sub-steps as the state transition of the exact model (see Section 6.2.2). However, the number of sub-steps reduces to three since only order processing at one process, namely the aggregated process, has to be modelled.

Due to these simplifications, state transition and transition probability can be expressed in a closed formula as follows

$$P(\mathbf{X}^{t+1} = \mathbf{z} \mid \mathbf{X}^t = \mathbf{x}) = \sum\_{(\mathbf{g}, \bar{h}) \in \mathcal{L}(\mathbf{x}, \mathbf{z})} P(\mathbf{G} = \mathbf{g}) \cdot P(\bar{H} = \bar{h})$$

$$\mathcal{L}(\mathbf{x}, \mathbf{z}) = \left\{ (\mathbf{g}, \bar{h}) \in \mathcal{G} \times \bar{\mathcal{H}} \mid \right.$$

$$g\_k + \max \left\{ 0; x\_{k+1} - \max \left\{ 0; \bar{h} - \sum\_{l=-R}^{k} x\_l \right\} \right\} = z\_k \qquad (7.6)$$

$$\forall k \in \mathcal{K} \; \middle\mid \left\{ e\_{max} \right\}$$

$$\wedge \; g\_{e\_{max}} = z\_{e\_{max}} \right\} \qquad \forall \mathbf{x}, \mathbf{z} \in \mathcal{X}.$$

The number of unprocessed orders z<sup>k</sup> with a due date of k time periods at the beginning of time period (t + 1) is the sum of the number of incoming orders per time period g<sup>k</sup> with a lead time of k time periods at the beginning of time period (t + 1) and the number of unprocessed orders with a due date of (k + 1) time periods at the end of time period t

$$\max\left\{0; x\_{k+1} - \max\left\{0; \bar{h} - \sum\_{l=-R}^{k} x\_l\right\}\right\}.$$

The number of unprocessed orders with a due date of (k + 1) time periods at the end of time period t is either zero or it corresponds to the difference of the number of unprocessed orders xk+1 with a due date of (k + 1) time periods at the beginning of time period t and the residual processing performance of time period t remaining after all orders with a due date of (l < k + 1) time periods have already been processed. The number of unprocessed orders zemax with a due date of emax time periods at the beginning of time period (t + 1) equals the number of incoming orders per time period gemax with a lead time of emax time periods at the beginning of time period (t + 1).

#### **7.2.3 State Space**

We derive the state space X of the Markov chain based on the same considerations as in the exact model (see Section 6.2.3). The state space X is non-negative, and its upper bound

$$\mathbf{O} = \begin{pmatrix} O\_{-R} & \dots & O\_{e\_{max}} \end{pmatrix} \tag{7.7}$$

is defined as follows

$$\begin{aligned} O\_k &= \left(e\_{\max} + 1\right) \cdot a\_{\max} & \forall k \in \{-R, \dots, -1\} \\ O\_k &= \left(e\_{\max} - k + 1\right) \cdot a\_{\max} & \forall k \in \{0, \dots, e\_{\max}\} \end{aligned} \tag{7.8}$$

#### **7.2.4 Limiting Distribution**

The Markov chain has a finite state space (see Section 7.2.3), and it is possible to reach every state of the Markov chain from every other state either by a singlestep transition or by an indirect transition via a finite number of other states. In the case of unreachable states, these states are excluded from the computations (see Section 10.1.1). We then find an irreducible subset of the state space. We assume that there is a tuple (g, h˜) ∈ G × H˜ for which the total number of incoming orders per time period equals the processing performance per time period:

$$\sum\_{k \in \mathcal{K}} g\_k = \bar{h}.$$

Thus, the irreducible subset of the state space contains state x with

$$x\_k = g\_k \quad \forall k \in \mathcal{K}$$

$$\sum\_{k \in \mathcal{K}} x\_k = \bar{h},$$

for which

$$P(\mathbf{X}^{t+1} = \mathbf{x} \mid \mathbf{X}^t = \mathbf{x}) = P(\bar{H} = \bar{h}) \cdot P(\mathbf{G} = \mathbf{g}) > 0$$

holds. Consequently, state x as well as all other states of the irreducible subset of the state space are aperiodic. In conclusion, the limiting distribution π˜ of the Markov chain exists, and we calculate the limiting distribution π˜ of a given initial state as described in Section 5.2.2.

## **7.3 Performance Measures of the Order Fulfilment System**

In the simplified model, we calculate the same performance measures of the order fulfilment system as in the exact model (see Table 6.1) using the formulas introduced in Section 6.3. However, the subsequent adjustments are necessary:


In the exact model, we subdivide the number of unprocessed orders into two key figures regarding their time of consideration (see Section 6.3.1). This differentiation is not meaningful in the simplified model since in single-stage order fulfilment systems, the point in time immediately before the start of order processing at a process corresponds to the beginning of the time period. Furthermore, the differentiation between process- and system-specific key figures regarding the number of unprocessed orders, the number of lost sales, and the number of processed orders is no longer reasonable in the simplified model.

## **7.4 Memory and Computation Time Requirements**

Memory and computation time requirements of the simplified model for performance analysis depend on the same factors as the ones of the exact model for performance analysis (see Section 6.4): The size of the state space of the Markov chain and the discretisation of the stochastic parameters.

In the simplified model, the size of the state space |X | of the Markov chain is derived from the lower and upper limits of the state space X (see Section 7.2.3) as follows

$$\begin{aligned} |\mathcal{X}| &= \prod\_{k \in \mathcal{K}} (O\_k + 1) \\ |\mathcal{X}| &= \left[ (e\_{\max} + 1) \cdot a\_{\max} + 1 \right]^R \cdot \prod\_{k=0}^{e\_{\max}} \left( (e\_{\max} - k + 1) \cdot a\_{\max} + 1 \right) . \end{aligned} \tag{7.9}$$

An upper limit of the size of state space is obtained by the following estimation

$$|\mathcal{X}| \le \left[ (e\_{\max} + 1) \cdot a\_{\max} + 1 \right]^{(R + c\_{\max} + 1)}.\tag{7.10}$$

From equation (7.10), we recognize that the maximum backlog duration R and the maximum possible values of the number of incoming order per time period A and the lead time E determine the size of the state space of the Markov chain. The maximum backlog duration R and the maximum possible lead time emax have a major impact on the size of the state space since they determine the exponent in equation (7.10).

## **7.5 Chapter Conclusion**

In this chapter, we introduced a simplified analytical model for performance analysis of multi-stage, stochastic order fulfilment systems with levelled order release and customer-required order deadlines by


The simplified model can be seen as a special case of the exact model introduced in Chapter 6 since it ensues from the exact model when modelling is restricted to single-stage order fulfilment systems. In Chapter 8, we will evaluate and compare both models for performance analysis regarding modelling accuracy, accuracy of performance analysis, and memory and computation time requirements.

# **8 Evaluation of Models for Performance Analysis**

The exact and the simplified model for performance analysis differ regarding the modelling of partially processed orders that are transmitted between the processing steps of a multi-stage order fulfilment system. The exact model provides an exact recording of partially processed orders, whereas partially processed orders are neglected in the simplified model. This chapter aims at evaluating and comparing the exact and the simplified model for performance analysis regarding modelling accuracy (see Section 8.1), accuracy of performance analysis (see Section 8.2), and memory and computation time requirements (see Section 8.3). By summarising the results of this chapter, Section 8.4 provides recommendations when to use which model for performance analysis of multi-stage, stochastic order fulfilment systems with levelled order release and customer-required deadlines.

## **8.1 Modelling Accuracy**

The main difference between the exact and the simplified model concerns the modelling of partially processed orders that are transmitted between the processing steps of a multi-stage order fulfilment system. The exact model considers partially processed orders by exactly recording the unprocessed orders per process (see Figure 5.1). In contrast, partially processed orders are neglected in the simplified model (see Figure 5.2) since an order is either completely unprocessed or completely processed. In the following, we initially derive two hypotheses on the impact of the modelling of partially processed orders on the modelling accuracy of multi-stage, stochastic order fulfilment systems with levelled order release and customer-required deadlines (see Section 8.1.1). Subsequently, we prove that these hypotheses hold (see Section 8.1.2).

#### **8.1.1 Hypotheses**

In order fulfilment systems with low utilisation, the expected value of the number of incoming orders per time period is significantly smaller than the expected value of the aggregated processing performance per time period (see equation (7.4)). The probability that the number of incoming orders per time period exceeds the minimum of the processing performances per time period of the processes p ∈ P in any time period is close to zero:

$$\begin{split} \frac{E(A)}{E(\bar{H})} &< <1 \Leftrightarrow E(A) <  \min\_{p \in \mathcal{P}} \{h\_p\} \right) \to 0. \end{split} (8.1)$$

Consequently, the majority of the incoming orders per time period is processed immediately after their time of arrival without being buffered at any processing step, and the probability that partially processed orders occur at any processing step is close to zero.

In contrast, in order fulfilment systems with high utilisation, the expected value of the number of incoming orders per time period is approximately as high as the expected value of the aggregated processing performance per time period. There is a significant probability that the number of incoming orders per time period exceeds the minimum of the processing performances per time period of the processes p ∈ P in any time period:

$$\begin{split} \frac{E(A)}{E(\bar{H})} &\to 1 \Leftrightarrow E(A) \approx E(\bar{H}) \\ &\Leftrightarrow \sum\_{a \in \mathcal{A}} a \cdot P(A = a) \approx \sum\_{\tilde{h} \in \tilde{\mathcal{H}}} \bar{h} \cdot P(\bar{H} = \bar{h}) \\ &\Leftrightarrow \sum\_{\mathbf{g} \in \mathcal{G}} \left( \sum\_{k \in \mathcal{K}} g\_k \right) P(\mathbf{G} = \mathbf{g}) \approx \sum\_{\mathbf{h} \in \mathcal{H}} \min\_{p \in \mathcal{P}} \{h\_p\} \cdot P(\mathbf{H} = \mathbf{h}) \\ &\Rightarrow \sum\_{(\mathbf{g}, \mathbf{h}) \in \mathcal{G} \times \mathcal{H}} P\left(\sum\_{k \in \mathcal{K}} g\_k > \min\_{p \in \mathcal{P}} \{h\_p\} \right) >> 0. \end{split} (8.2)$$

Consequently, there is a significant probability that temporary buffers occur at any processing step of the multi-stage system. Buffers in a multi-stage system contain partially processed orders. In conclusion, we state that partially processed orders predominantly occur in systems with high utilisation. Hence, the focus for investigating the impact of the modelling of partially processed orders is on order fulfilment systems with high utilisation in the following.

In order fulfilment systems with high utilisation, the processing performance per time period of the system is the limiting factor of system throughput (see equation (8.2)). Regarding the simplified model, this means that system throughput depends on the aggregated processing performance per time period, and thus on the minimum of the processing performances per time period of processes p ∈ P (see equation (7.1)). In the exact model, the relationship is more complicated since there is no parameter describing the processing performance per time period of the whole system. Instead, system throughput of systems with high utilisation depends on


To illustrate these dependencies, we consider the processing step pmax of a multi-stage system in any time period in detail, without loss of generality: If the processing performance per time period hpmax of process pmax is at most as high as the processing performance per time period h(pmax−1) of process(pmax − 1), (hpmax ≤ hpmax−1), the number of completed orders equals the processing performance per time period hpmax of process pmax independent of the number of partially processed orders at process pmax. Otherwise, if the processing performance per time period hpmax of process pmax exceeds the processing performance per time period h(pmax−1) of process(pmax − 1), (hpmax > hpmax−1), the number of completed orders equals the processing performance per time period h(pmax−1) of process (pmax − 1) if no partially processed orders are buffered at process pmax. However, it depends on the processing performance per time period hpmax of process pmax if there is a sufficiently high number of partially processed orders buffered at process pmax. The underlying precondition is that temporary buffers of partially processed orders can occur at any processing step of the system. Referring to common bottleneck definitions in literature (e.g. Lawrence and Buss (1994), Lai et al. (2021)), we call them systems with a shifting bottleneck.

A system with a finite set of processes P and a processing performance per time period H is said to be a system with a shifting bottleneck if the processing performance per time period H meets the following conditions:

$$\exists \text{ } \mathbf{h} \in \mathcal{H}: h\_i > h\_j \qquad i, j \in \mathcal{P}, i \neq j \tag{8.3}$$

$$\exists \text{ } \mathbf{h} \in \mathcal{H}: h\_i \le h\_j \qquad i, j \in \mathcal{P}, i \ne j. \tag{8.4}$$

In contrast, it is said to be a system with a static bottleneck if there is a process p <sup>∗</sup> ∈ P, called static bottleneck, for which the following condition holds:

$$\forall \mathbf{h} \in \mathcal{H}: h\_{p^\*} < h\_j \qquad j \in \mathcal{P}, j \neq p^\*. \tag{8.5}$$

In the special case of systems with a static bottleneck, buffers of partially processed orders only occur at the bottleneck and at processes upstream the bottleneck. In this case, the number of completed orders per time period does not depend on the three factors mentioned above, but it only depends on the processing performance per time period of the bottleneck. However, the processing performance per time period of the bottleneck equals the minimum of the processing performances per time period of processes p ∈ P (see equation (8.5)), and thus we observe the same relation as in the simplified model. In conclusion, we state the following: In systems with high utilisation and a static bottleneck, system throughput only depends on the processing performance per time period of the bottleneck, and modelling of partially processed orders can be neglected. Otherwise, in systems with high utilisation and a shifting bottleneck, system throughput additionally depends on the temporary number and location of partially processed orders. These findings lead to the following hypothesis:

Hypothesis 1. For performance analysis of order fulfilment systems with high utilisation and a shifting bottleneck, there is a loss of modelling accuracy when using the simplified model instead of the exact model.

As mentioned before, in the simplified model, system throughput of systems with high utilisation is determined by the aggregated processing performance per time period, which corresponds to the minimum of the processing performances per time period of processes p ∈ P (see equation (7.1)). In contrast, in the exact model, system throughput of systems with high utilisation and a shifting bottleneck is determined by the processing performance per time period of process pmax, when assuming a sufficiently high number of partially processed orders buffered at process pmax. There are realisations h ∈ H for which the processing performance per time period of process pmax is the minimum of the processing performances per time period of processes p ∈ P (see equation (8.4)), but there are also realisations h ∈ H for which this condition does not hold (see equation (8.3)). These findings lead to the following hypothesis:

Hypothesis 2. In systems with high utilisation and a shifting bottleneck, the simplified model underestimates the exact value of the expected value of the system throughput.

### **8.1.2 Mathematical Proof**

In the following, we prove hypotheses 1 and 2. For this, some preliminary remarks on the characteristics of order processing at any process p ∈ P (see Section 8.1.2.1) and the lower bound of the total number of unprocessed orders at process p ∈ P immediately before the start of order processing at process p (see Section 8.1.2.2) are required. Subsequently, we derive a formula of system throughput for each system configuration – system with a static and a shifting bottleneck – in each modelling approach – exact and simplified model – under the precondition of high system utilisation (see equation (8.2)) (see Sections 8.1.2.3- 8.1.2.5). By comparing these formulas in Section 8.1.2.6, we finally prove that hypotheses 1 and 2 hold.

#### **8.1.2.1 Preliminaries on Order Processing at Process** p ∈ P

Order processing at process p in time period t depends on the realisation H<sup>p</sup> = h<sup>p</sup> of the processing performance per time period of process p in time period t. Following the principles of levelled order release, order processing at process p starts with orders having a due date of (−R) time periods. Thus, in time period t, a processing performance of h<sup>p</sup> orders is available for processing orders with a due date of (−R) time periods. In contrast, for processing orders with a due date of(k > −R)time periods, only the so-called residual processing performance h res p is available that corresponds to the processing performance remaining after all orders with a due date of (l < k) time periods have already been processed.

The residual processing performance h res <sup>p</sup> of process p in time period t can be written as a function of due date k ∈ K

$$h\_p^{res}(k) = \max\left\{0; h\_p - \sum\_{l=-R}^{k-1} y\_{p,l}^{(p-1)}\right\},\tag{8.6}$$

.

whereby h res p (k) specifies the processing performance remaining for processing orders with a due date of k time periods at process p in time period t. Due to the principles of levelled order release, h res p (k) is monotonously decreasing in k:

$$h\_p^{res}(-R) = \max\left\{0; h\_p\right\} = h\_p$$

$$h\_p^{res}(-R+1) = \max\left\{0; h\_p - y\_{p,-R}^{(p-1)}\right\}$$

$$\vdots$$

$$h\_p^{res}(e\_{max}) = \max\left\{0; h\_p - \sum\_{l=-R}^{e\_{max}-1} y\_{p,l}^{(p-1)}\right\}$$

We define the critical due date k ∗ <sup>p</sup> of process p in time period t as the shortest due date whose orders are not processed in time period t:

$$\begin{aligned} k\_p^\* &= \begin{cases} \min\left\{ k \in \mathcal{K} \mid h\_p^{res}(k) = 0 \right\} & \exists k \in \mathcal{K} : h\_p^{res}(k) = 0\\ \infty & \text{otherwise} \end{cases} \\ k\_p^\* &= \begin{cases} \min\left\{ k \in \mathcal{K} \mid \sum\_{l=-R}^{k-1} y\_{p,l}^{(p-1)} \ge h\_p \right\} & \exists k \in \mathcal{K} : h\_p^{res}(k) = 0\\ \infty & \text{otherwise} \end{cases} \end{aligned} \tag{8.7}$$

The critical due date k ∗ <sup>p</sup> of process p in time period t can be interpreted as follows: No order with a due date of (k ≥ k ∗ p ) time periods is processed at process p in time period t. In contrast, unprocessed orders with a due date of (k ∗ <sup>p</sup> − 1) time periods are either partially or completely processed at process p in time period t, and all unprocessed orders having a due date of (k ≤ k ∗ <sup>p</sup> − 2) time periods are completely processed at process p in time period t.

The number of processed orders per time period mp,k at process p with a due date of k time periods at their time of processing in time period t depends on number of unprocessed orders at process p with a due date of k time periods immediately before the start of order processing at process p and the residual processing performance (see equation (6.31)):

$$\begin{split} m\_{p,k} &= \min \left\{ y\_{p,k}^{(p-1)}; \max \left\{ 0; h\_p - \sum\_{l=-R}^{k-1} y\_{p,l}^{(p-1)} \right\} \right\} \\ m\_{p,k} &\overset{(8.6)}{=} \min \left\{ y\_{p,k}^{(p-1)}; h\_p^{res}(k) \right\} \\ m\_{p,k} &\overset{(8.7)}{=} \begin{cases} y\_{p,k}^{(p-1)} & k \in \mathcal{K}, k \le k\_p^\* - 2 \\ h\_p - \sum\_{l=-R}^{k-1} y\_{p,l}^{(p-1)} & k \in \mathcal{K}, k = k\_p^\* - 1 \\ 0 & k \in \mathcal{K}, k \ge k\_p^\*. \end{cases} \end{split} \tag{8.8}$$

Consequently, the total number of processed orders per time period f<sup>p</sup> at process p in time period t (see equation (6.33)) can be written as follows

$$\begin{split} f\_p &= \sum\_{k \in \mathcal{K}} m\_{p,k} \\ f\_p &\overset{(8.8)}{=} \sum\_{k=-R}^{k\_p^\*-2} y\_{p,k}^{\{p-1\}} + \left( h\_p - \sum\_{l=-R}^{k\_p^\*-2} y\_{p,l}^{\{p-1\}} \right) + \sum\_{k=k\_p^\*}^{e\_{max}} 0 \\ f\_p &\overset{(8.7)}{=} \begin{cases} \sum\_{k \in \mathcal{K}} y\_{p,k}^{\{p-1\}} & k\_p^\* = \infty \\ h\_p & k\_p^\* \in \mathcal{K}. \end{cases} \end{split} \tag{8.9}$$

#### **8.1.2.2 Preliminaries on Lower Bound of Total Number of Unprocessed Orders at Process** p ∈ P **Immediately Before the Start of Order Processing**

Using the notation of the state transition of the Markov chain (see Section 6.2.2), y (p−1) <sup>p</sup> , p ∈ P, denotes the number of unprocessed orders at process p ∈ P immediately before the start of order processing at process p in any time period t.

The total number of unprocessed orders P <sup>k</sup>∈K y (0) 1,k at process p = 1 immediately before the start of order processing at this process corresponds to the total number of unprocessed orders at process p = 1 at the beginning of time period t. The last sub-step of every state transition specifies the order income per time period G = g of the following time period (see Section 6.2.2). Thus, the total number of unprocessed orders at process p = 1 immediately before the start of order processing at this process is at least as high as the total number of incoming orders at the beginning of time period t:

$$\sum\_{k \in \mathcal{K}} y\_{1,k}^{(0)} \ge \sum\_{k \in \mathcal{K}} g\_k. \tag{8.10}$$

The total number of unprocessed orders P <sup>k</sup>∈K y (p−1) p,k at process p ∈ P \ {1} immediately before the start of order processing at this process in time period t is the sum of the number of unprocessed orders at process p at the end of the previous time period (t − 1) and the number of processed orders at the previous processing step (p − 1) in time period t (see equation (6.10))

$$\begin{split} \sum\_{k \in \mathcal{K}} y\_{p,k}^{(p-1)} &= \sum\_{k \in \mathcal{K}} y\_{p,k}^{(p-2)} \\ &+ \sum\_{k \in \mathcal{K}} \min \left\{ y\_{(p-1),k}^{(p-2)}; \max \left\{ 0; h\_{p-1} - \sum\_{l=-R}^{k-1} y\_{(p-1),l}^{(p-2)} \right\} \right\} \\ \sum\_{k \in \mathcal{K}} y\_{p,k}^{(p-1)} &\overset{(8.9)}{=} \begin{cases} \sum\_{k \in \mathcal{K}} y\_{p,k}^{(p-2)} + \sum\_{k \in \mathcal{K}} y\_{(p-1),k}^{(p-2)} & k\_{p-1}^\* = \infty \\ \sum\_{k \in \mathcal{K}} y\_{p,k}^{(p-2)} + h\_{p-1} & k\_{p-1}^\* \in \mathcal{K}, \end{cases} \end{split} \tag{8.11}$$

whereby H = h specifies the processing performance per time period in time period t. Based on equations (8.10) and (8.11), we can derive a lower bound for the total number of unprocessed orders at any process p ∈ P \ {1} immediately before the start of order processing at this process in time period t. For process p = 2, the lower bound is derived as follows

$$\begin{split} \sum\_{k \in \mathcal{K}} y\_{2,k}^{(1)} &= \begin{cases} \sum\_{k \in \mathcal{K}} y\_{2,k}^{(0)} + \sum\_{k \in \mathcal{K}} y\_{1,k}^{(0)} & k\_1^\* = \infty \\ \sum\_{k \in \mathcal{K}} y\_{2,k}^{(0)} + h\_1 & k\_1^\* \in \mathcal{K} \end{cases} \\ \overset{(8.10)}{\Rightarrow} \sum\_{k \in \mathcal{K}} y\_{2,k}^{(1)} &\geq \begin{cases} \sum\_{k \in \mathcal{K}} g\_k & k\_1^\* = \infty \\ h\_1 & k\_1^\* \in \mathcal{K}. \end{cases} \end{split} \tag{8.12}$$

For process p = 3, the lower bound is derived as follows

$$\sum\_{k \in \mathcal{K}} y\_{3,k}^{(2)} = \begin{cases} \sum\_{k \in \mathcal{K}} y\_{3,k}^{(1)} + \sum\_{k \in \mathcal{K}} y\_{2,k}^{(1)} & k\_2^\* = \infty \\ \sum\_{k \in \mathcal{K}} y\_{3,k}^{(1)} + h\_2 & k\_2^\* \in \mathcal{K} \end{cases}$$

$$y\_{3,k}^{(8.12)} \sum\_{k \in \mathcal{K}} y\_{3,k}^{(2)} \ge \begin{cases} \sum\_{k \in \mathcal{K}} g\_k & k\_1^\* = \infty, k\_2^\* = \infty \\ h\_1 & k\_1^\* \in \mathcal{K}, k\_2^\* = \infty \\ h\_2 & k\_2^\* \in \mathcal{K}. \end{cases} \tag{8.13}$$

Consequently, the lower bound of any process p ∈ P \ {1} is defined as follows

$$\sum\_{k \in \mathcal{K}} y\_{p,k}^{(p-1)} \ge \begin{cases} \sum\_{k \in \mathcal{K}} g\_k & k\_1^\*, \dots, k\_{p-1}^\* = \infty \\ h\_j & k\_j^\* \in \mathcal{K}, k\_{j+1}^\*, \dots, k\_{p-1}^\* = \infty, \\ & j \in \{1, \dots, (p-1)\}. \end{cases} \tag{8.14}$$

#### **8.1.2.3 System Throughput in Systems With a Static Bottleneck in the Exact Model**

Under the precondition of high utilisation (see equation (8.2)), in systems with a static bottleneck p ∗ (see equation (8.5)), the lower bound of the total number of unprocessed orders at the bottleneck p ∗ immediately before the start of order processing at process p ∗ (see equation (8.14)) simplifies to

$$\sum\_{k \in \mathcal{K}} y\_{p^\*,k}^{(p^\*-1)} > h\_{p^\*} \qquad \forall \mathbf{h} \in \mathcal{H}.\tag{8.15}$$

Based on this, we conclude that the critical due date k ∗ <sup>p</sup><sup>∗</sup> of the bottleneck p ∗ satisfies k ∗ <sup>p</sup><sup>∗</sup> ∈ K and that the throughput fp<sup>∗</sup> of the bottleneck p ∗ is given by its processing performance per time period hp<sup>∗</sup> for all realisations h ∈ H:

$$\sum\_{k \in \mathcal{K}} y\_{p^{\bullet}, k}^{(p^{\bullet}-1)} > h\_{p^{\bullet}} \quad \forall \mathbf{h} \in \mathcal{H} \stackrel{(8.7)}{\Rightarrow} k\_{p^{\bullet}}^{\bullet} \in \mathcal{K} \qquad \forall \mathbf{h} \in \mathcal{H} \tag{8.16}$$

$$\stackrel{(8.9)}{\Rightarrow} f\_{p^{\bullet}} = h\_{p^{\bullet}} \quad \forall \mathbf{h} \in \mathcal{H}.$$

If the last processing step pmax is the bottleneck of the system, system throughput directly results from equation (8.16) as follows

$$f \stackrel{(6.35)}{=} f\_{p\_{max}} = f\_{p^\*} \stackrel{(8.16)}{=} h\_{p^\*} \qquad \forall \mathbf{h} \in \mathcal{H}. \tag{8.17}$$

Otherwise, if any other processing step p ∈ P \ {pmax} represents the bottleneck p <sup>∗</sup> of the order fulfilment system, we additionally consider order processing at the subsequent processing step (p <sup>∗</sup> + 1) of the bottleneck p ∗ . Since the throughput of the bottleneck p ∗ corresponds to its processing performance per time period hp<sup>∗</sup> for all realisations h ∈ H (see equation (8.16)), and the processing performance per time period h(p∗+1) of process (p <sup>∗</sup> + 1) is higher than the one of the bottleneck p ∗ for all realisations h ∈ H (see equation (8.5)), there never remain any unprocessed order at process (p <sup>∗</sup> + 1) at the end of any time period. The total number of unprocessed orders at process (p <sup>∗</sup> + 1) immediately before the start of order processing at process(p <sup>∗</sup> + 1)(see equation (8.11)) simplifies to

$$\sum\_{k \in \mathcal{K}} y\_{(p^\star + 1), k}^{(p^\star)} = h\_{p^\star} \qquad \forall \mathbf{h} \in \mathcal{H}.\tag{8.18}$$

Based on this, we conclude that the throughput f(p∗+1) of process (p <sup>∗</sup> + 1) is given by the processing performance per time period hp<sup>∗</sup> of the bottleneck p ∗ for all realisations h ∈ H:

$$\sum\_{k \in \mathcal{K}} y\_{(p^\* + 1), k}^{(p^\*)} = h\_{p^\*} \quad \forall \mathbf{h} \in \mathcal{H} \stackrel{(8.5)}{\Rightarrow} \sum\_{k \in \mathcal{K}} y\_{(p^\* + 1), k}^{(p^\*)} < h\_{(p^\* + 1)} \quad \forall \mathbf{h} \in \mathcal{H}$$

$$\stackrel{(8.7)}{\Rightarrow} k\_{(p^\* + 1)}^{\*} = \infty \qquad \forall \mathbf{h} \in \mathcal{H} \tag{8.19}$$

$$\stackrel{(8.9)}{\Rightarrow} f\_{(p^\* + 1)} = \sum\_{k \in \mathcal{K}} y\_{(p^\* + 1), k}^{(p^\*)} \quad \forall \mathbf{h} \in \mathcal{H}$$

$$\stackrel{(8.18)}{\Rightarrow} f\_{(p^\* + 1)} = h\_{p^\*} \quad \forall \mathbf{h} \in \mathcal{H}.$$

The same reasoning holds for all subsequent processes p ∈ {(p <sup>∗</sup> + 2), . . . , pmax}, and thus system throughput is given by the processing performance per time period hp<sup>∗</sup> of the bottleneck p ∗

$$f = h\_{p^\*} \qquad \forall \mathbf{h} \in \mathcal{H}.\tag{8.20}$$

In conclusion, we showed that under the precondition of high utilisation (see equation (8.2)), system throughput f of order fulfilment systems with a static bottleneck p ∗ in any time period is determined by the processing performance per time period hp<sup>∗</sup> of the bottleneck p ∗ for all realisations h ∈ H, independent of the position of the bottleneck p <sup>∗</sup> ∈ P

$$f = h\_{p^\*} \stackrel{(8.5)}{=} \min\_{p \in \mathcal{P}} \{h\_p\} \qquad \forall \mathbf{h} \in \mathcal{H}. \tag{8.21}$$

#### **8.1.2.4 System Throughput in Systems With a Shifting Bottleneck in the Exact Model**

For the analysis of system throughput in systems with a shifting bottleneck, we focus on the last processing step pmax of the multi-stage system in any time period t, whereby H = h specifies the processing performance per time period in time period t. The following considerations assume high utilisation (see equation (8.2)).

If the last processing step is the bottleneck of the system in time period t (see equation (8.4)), following the reasoning of equation (8.16), system throughput f in time period t corresponds to the minimum of the processing performances per time period of processes p ∈ P (see equation (8.17)).

Otherwise, if the last processing step is not the bottleneck of the system in time period t (see equation (8.3)), the throughput f<sup>p</sup>max of process pmax in time period t either corresponds to the processing performance per time period of process pmax in time period t or to the total number of unprocessed orders at process pmax immediately before the start of order processing at process pmax in time period t, depending on the value of the critical due date k ∗ pmax of process pmax in time period t (see equation (8.9)):

• Critical due date of k ∗ pmax ∈ K

$$f\_{p\_{max}} = h\_{p\_{max}} \stackrel{(8.3)}{\Rightarrow} f\_{p\_{max}} > \min\_{p \in \mathcal{P}} \{h\_P\}.\tag{8.22}$$

• Critical due date of k ∗ <sup>p</sup>max = ∞

$$\begin{split} f\_{p\_{\max}} &= \sum\_{k \in \mathcal{K}} y\_{p\_{\max}, k}^{(p\_{\max}-1)} \\ \stackrel{(8.14)}{\Rightarrow} f\_{p\_{\max}} &\geq \begin{cases} \sum\_{k \in \mathcal{K}} g\_k & k\_1^\*, \dots, k\_{p\_{\max}-1}^\* = \infty \\ h\_j & k\_j^\* \in \mathcal{K}, k\_{j+1}^\*, \dots, k\_{p\_{\max}-1}^\* = \infty, \\ &\qquad j \in \{1, \dots, (p\_{\max}-1)\} \end{cases} \\ \stackrel{(8.3)}{\Rightarrow} f\_{p\_{\max}} &\geq \min\_{p \in \mathcal{P}} \{h\_p\}. \end{split} \tag{8.23}$$

In conclusion, we showed that under the precondition of high utilisation (see equation (8.2)), system throughput f of order fulfilment systems with a shifting bottleneck in any time period depends on the realisation h ∈ H as follows

$$f = \min\_{p \in \mathcal{P}} \{ h\_p \} \qquad \forall \mathbf{h} \in \mathcal{H} : h\_{p\_{max}} = \min\_{p \in \mathcal{P}} \{ h\_p \} \tag{8.24}$$

$$f > \min\_{p \in \mathcal{P}} \{ h\_p \} \qquad \forall \mathbf{h} \in \mathcal{H}: h\_{p\_{max}} > \min\_{p \in \mathcal{P}} \{ h\_p \}, k\_{p\_{max}}^\* \in \mathcal{K} \tag{8.25}$$

$$f \ge \min\_{p \in \mathcal{P}} \{ h\_p \} \qquad \forall \mathbf{h} \in \mathcal{H} : h\_{p\_{max}} > \min\_{p \in \mathcal{P}} \{ h\_p \}, \\ k\_{p\_{max}}^\* = \infty. \tag{8.26}$$

#### **8.1.2.5 System Throughput in the Simplified Model**

In the simplified model, system throughput in time period t either corresponds to the aggregated processing performance per time period h˜ in time period t or to the total number of unprocessed orders at the beginning of time period t, depending on the value of the critical due date k ∗ in time period t (simplification of equation (8.9)):

$$f = \begin{cases} \bar{h} & k^\* \in \mathcal{K} \\ \sum\_{k \in \mathcal{K}} y\_k^{(0)} & k^\* = \infty. \end{cases} \tag{8.27}$$

Under the precondition of high utilisation (see equation (8.2)), we conclude from the lower bound of the total number of unprocessed orders (see equation (8.10)) that system throughput f in any time period is determined by the aggregated processing performance per time period h˜ for all realisations h˜ ∈ H˜:

$$\sum\_{k \in \mathcal{K}} y\_k^{(0)} \ge \sum\_{k \in \mathcal{K}} g\_k \stackrel{(8.2)}{\Rightarrow} \sum\_{k \in \mathcal{K}} y\_k^{(0)} > \bar{h} \qquad \forall \bar{h} \in \bar{\mathcal{H}}$$

$$\stackrel{(8.7)}{\Rightarrow} k^\* \in \mathcal{K} \qquad \forall \bar{h} \in \bar{\mathcal{H}} \tag{8.28}$$

$$\stackrel{(8.27)}{\Rightarrow} f = \bar{h} \qquad \forall \bar{h} \in \bar{\mathcal{H}}$$

$$\stackrel{(7.1)}{\Rightarrow} f = \min\_{p \in \mathcal{P}} \{h\_p\} \quad \forall \mathbf{h} \in \mathcal{H}.$$

#### **8.1.2.6 Comparison of System Throughput in the Exact and the Simplified Model**

Table 8.1 summarises the values of system throughput in time period t for each system configuration – system with a static or a shifting bottleneck – in each modelling approach – exact or simplified model – under the precondition of high system utilisation (see equation (8.2)) as derived in the previous sections.


Table 8.1: Comparison of system throughput f in any time period t in the exact and the simplified model (under the precondition of high system utilisation).

We proved that system throughput of order fulfilment systems with a static bottleneck in any time period t equals the minimum of the processing performances per time period of processes p ∈ P under the precondition of high utilisation. It is calculated exactly based on the simplified model. There is no loss of modelling accuracy for performance analysis of order fulfilment systems with a static bottleneck when using the simplified model, which does not model partially processed orders, instead of the exact model. Consequently, the modelling of partially processed orders has no impact on the modelling accuracy of order fulfilment systems with a static bottleneck.

Furthermore, we proved that system throughput of order fulfilment systems with a shifting bottleneck in any time period t is at least as high as the minimum of the processing performances per time period of processes p ∈ P under the precondition of high utilisation. For systems with high utilisation and a shifting bottleneck, the simplified model only provides a lower bound of the exact system throughput since it approximates system throughput in any time period t by the minimum of the processing performances per time period of processes p ∈ P. There are modelling inaccuracies for performance analysis of order fulfilment systems with high utilisation and a shifting bottleneck when using the simplified model instead of the exact one. Consequently, the modelling of partially processed orders affects modelling accuracy of order fulfilment systems with high utilisation and a shifting bottleneck. The expected value of the system throughput calculated based on the simplified model is smaller than the one calculated based on the exact model.

By these findings, hypotheses 1 and 2 are mathematically proven.

## **8.2 Accuracy of Performance Analysis**

For performance analysis of an order fulfilment system, several performance measures (see Table 6.1) are derived from the limiting distribution of the Markov chain in both the exact and the simplified model. Due to the proven modelling inaccuracies of the simplified model (see Section 8.1), we expect inaccuracies of performance analysis when using the simplified model. In the following, we initially derive two hypotheses on the impact of the modelling inaccuracies of the simplified model on the accuracy of performance analysis (see Section 8.2.1). Subsequently, we verify these hypotheses in a numerical analysis (see Section 8.2.2).

### **8.2.1 Hypotheses**

In Section 8.1, we showed that the simplified model suffers from modelling inaccuracies for performance analysis of order fulfilment systems with high utilisation and a shifting bottleneck indicated by deviations of system throughput in any time period t. Since the customer-related performance measures of the order fulfilment system are derived from the system throughput and the characteristics of completed orders, the results of Section 8.1 lead to the following hypothesis:

Hypothesis 3. For performance analysis of order fulfilment systems with high utilisation and a shifting bottleneck, the values of customer-related performance measures calculated based on the simplified model deviate from the exact ones calculated based on the exact model.

Furthermore, we proved in Section 8.1 that in systems with high utilisation and a shifting bottleneck, the expected value of the system throughput calculated based on the simplified model is smaller than the exact one calculated based on the exact model. Thus, fewer orders are completed on time, and the proportion of on-time completed orders on the total number of orders in the simplified model is smaller than the one in the exact model. Based on the definition of the service level (see Section 6.3.7), these findings lead to the following hypothesis:

Hypothesis 4. For performance analysis of order fulfilment systems with high utilisation and a shifting bottleneck, the simplified model underestimates the exact values of α-, β-, and γ-service level.

### **8.2.2 Numerical Results**

To verify hypotheses 3 and 4, we conduct a numerical analysis based on samples C and E. Sample C contains 256 two-stage order fulfilment systems, and sample E consists of 353 three-stage order fulfilment systems. The samples differ regarding the ranges of the parameters, especially regarding the discretisation of the number of incoming orders per time period A: Its range is given by A<sup>C</sup> = {3, 4, . . . , 14} in sample C and A<sup>E</sup> = {2, 3, . . . , 7} in sample E. A comprehensive description of the samples is given in Appendix C.2. Sample C consists of 167 systems with a shifting bottleneck (65% of sample size) and 89 systems with a static bottleneck. Sample E consists of 334 systems with a shifting bottleneck (95% of sample size) and 19 systems with a static bottleneck.

The numerical analysis focuses on α- and β-service level for simplicity reasons. As evaluation criteria, we calculate the absolute and the relative deviation of the values of α- and β-service level calculated based on the simplified model from the ones calculated based on the exact model, respectively.

#### **8.2.2.1 Order Fulfilment Systems With a Static Bottleneck**

Tables 8.2 and 8.3 present the absolute and the relative deviations of α- and β-service level between the simplified and the exact model for order fulfilment systems with a static bottleneck in samples C and E, respectively. Each sample is subdivided into multiple classes according to the utilisation of the order fulfilment system (see equation (6.7)). The tables give the minimum, the maximum, and the average values of the absolute and relative deviations for each class of utilisation as well as for the whole sample.

The absolute deviation of α-service level is between -5.17E-05 and 8.53E-05 in sample C and between -1.92E-05 and 2.56E-05 in sample E. We observe absolute deviations of β-service level between -3.14E-05 and 4.88E-05 in sample C and between -9.66E-06 and 2.05E-05 in sample E. Thus, the absolute values of the absolute deviations of α- and β-service level are smaller than 8.6E-05 in both


Table 8.2: Deviation of α-service level between the simplified and the exact model for order fulfilment systems with a static bottleneck.

<sup>1</sup> Difference between the value of α-service level calculated based on the simplified model and the one calculated based on the exact model.

<sup>2</sup> Difference between the value of α-service level calculated based on the simplified model and the one calculated based on the exact model, divided by the value calculated based on the exact model.

samples. We state these deviations to be negligible. They result from numerical errors during the calculations. By comparing the magnitude of the absolute and relative deviations of α- and β-service level of different classes of utilisation in each sample, we note that the utilisation of the order fulfilment system does not have any systematic impact on the absolute and relative deviations of neither αnor β-service level.

In conclusion, these results indicate that the chosen model for performance analysis – exact or simplified model – has a negligible impact on the values of α- and β-service level in order fulfilment systems with a static bottleneck.


Table 8.3: Deviation of β-service level between the simplified and the exact model for order fulfilment systems with a static bottleneck.

<sup>1</sup> Difference between the value of β-service level calculated based on the simplified model and the one calculated based on the exact model.

<sup>2</sup> Difference between the value of β-service level calculated based on the simplified model and the one calculated based on the exact model, divided by the value calculated based on the exact model.

#### **8.2.2.2 Order Fulfilment Systems With a Shifting Bottleneck**

Tables 8.4 and 8.5 present the absolute and the relative deviations of α- and β-service level between the simplified and the exact model for order fulfilment systems with a shifting bottleneck in samples C and E, respectively. Each sample is subdivided into multiple classes according to the utilisation of the order fulfilment system (see equation (6.7)). The tables give the minimum, the maximum, and the average values of the absolute and relative deviations for each class of utilisation as well as for the whole sample.


Table 8.4: Deviation of α-service level between the simplified and the exact model for order fulfilment systems with a shifting bottleneck.

<sup>1</sup> Difference between the value of α-service level calculated based on the simplified model and the one calculated based on the exact model.

<sup>2</sup> Difference between the value of α-service level calculated based on the simplified model and the one calculated based on the exact model, divided by the value calculated based on the exact model.

We observe absolute deviations of α-service level between -3.12E-02 and 7.85E-05 in sample C and between -2.78E-01 and 6.37E-05 in sample E. The absolute deviation of β-service level is between -6.69E-03 and 5.99E-05 in sample C and between -2.29E-01 and 5.43E-05 in sample E. These results indicate that for order fulfilment systems with a shifting bottleneck, there are data points in both samples for which the values of α- and β-service level calculated based on the simplified model deviate significantly from the ones calculated based on the exact model.


Table 8.5: Deviation of β-service level between the simplified and the exact model for order fulfilment systems with a shifting bottleneck.

<sup>1</sup> Difference between the value of β-service level calculated based on the simplified model and the one calculated based on the exact model.

<sup>2</sup> Difference between the value of β-service level calculated based on the simplified model and the one calculated based on the exact model, divided by the value calculated based on the exact model.

Figure 8.1 presents the relative deviations of α- and β-service level between the simplified and the exact model depending on the utilisation of the order fulfilment system for order fulfilment systems with a shifting bottleneck in samples C and E. It indicates that the utilisation of the order fulfilment system has a systematic impact on the relative deviations of α- and β-service level since the absolute value of the relative deviations of α- and β-service level increases as the utilisation converges towards one.

Figure 8.1: Deviation of α- and β-service level between the simplified and the exact model for order fulfilment systems with a shifting bottleneck.

In sample C, for systems with a utilisation smaller than 0.8, the relative deviation is between -0.003% and 0.005% for α-service level and between -0.001% and 0.004% for β-service level. In sample E, for systems with a utilisation smaller than 0.6, we observe a relative deviation of α-service level between -0.001% and 0.006% and of β-service level between -0.001% and 0.005%. Thus, for order fulfilment systems with a utilisation smaller than approximately 0.6, the absolute value of the relative deviation of α- and β-service level is at most 0.005% in both samples. We state these deviations to be negligible and conclude that the chosen model for performance analysis – exact or simplified model – has a negligible impact on the values of α- and β-service level in order fulfilment systems with a utilisation smaller than approximately 0.6 and with a shifting bottleneck.

In contrast, in sample C, the maximum absolute value of the relative deviation is 3.3% for α-service level and 1% for β-service level for systems with a utilisation of [0.9,1.0). In sample E, the maximum absolute value of the relative deviation is 17.3% for α-service level and 10.4% for β-service level when considering systems with a utilisation of [0.8,0.9). Furthermore, for systems with a utilisation of [0.9,1.0), we observe a maximum absolute value of the relative deviation of 47.6% for α-service level and 26.3% for β-service level. These results indicate that the values of α- and β-service level calculated based on the simplified model can strongly deviate from the ones calculated based on the exact model for systems with a utilisation higher than 0.8 and a shifting bottleneck.

By comparing the magnitude of the absolute and relative deviations of α- and β-service level between sample C and sample E, we observe that the number of processes of the considered order fulfilment system has a systematic impact on the deviations of α- and β-service level between the simplified and the exact model. The maximum absolute value of the relative deviation of α-service level is 3.3% for the two-stage systems of sample C, whereas it is 47.6% for the threestage systems of sample E. The maximum absolute value of the relative deviation of β-service level is 1% for the two-stage systems of sample C and 26.3% for the three-stage systems of sample E. These results indicate that the absolute values of the deviations of α- and β-service level between the simplified and the exact model increase with an increasing number of processes in the system.

### **8.2.3 Discussion**

For order fulfilment systems with a utilisation smaller than approximately 0.6 or systems with a static bottleneck, the chosen model for performance analysis – exact or simplified model – has a negligible impact on the values of αand β-service level in both samples. For these systems, the observed absolute deviations have a magnitude of at most 6.4E-05. In contrast, for order fulfilment systems with a utilisation higher than 0.8 and a shifting bottleneck, there are significant deviations of α- and β-service level between the simplified and the exact model. For instance, for systems with a utilisation of [0.9,1.0), we observe an absolute value of the relative deviation of α-service level of at most 3.3% and on average 0.5% in sample C and of at most 47.6% and on average 6.9% in sample E. In conclusion, these results confirm that the values of customer-related performance measures calculated based on the simplified model deviate from the exact ones calculated based on the exact model for performance analysis of order fulfilment systems with high utilisation and a shifting bottleneck (see hypothesis 3).

The deviations of α- and β-service level observed for order fulfilment systems with a utilisation higher than 0.8 and a shifting bottleneck are negative. Thus, the values of α- and β-service level calculated based on the simplified model are smaller than the ones calculated based on the exact model. These results confirm that the simplified model underestimates the exact values of α- and β-service level for performance analysis of order fulfilment systems with high utilisation and a shifting bottleneck (see hypothesis 4).

Furthermore, for order fulfilment systems with a utilisation higher than 0.6 and a shifting bottleneck, the number of processes of the considered order fulfilment system has a positive impact on the absolute magnitude of the deviations of αand β-service level between the simplified and the exact model. The modelling inaccuracies of the simplified model result from the fact that partially processed orders are not modelled in the simplified model (see Section 8.1). Partially processed orders occur at the buffers of the processes p ∈ P \ {1} of an order fulfilment system. Consequently, the higher the number of processes, the more partially processed orders can occur in the order fulfilment system, and the higher are the modelling inaccuracies and thus the inaccuracies of performance analysis of the simplified model.

In the context of workload balancing in order fulfilment systems, the focus is neither on systems with low utilisation nor on systems with a static bottleneck: On the one hand, in order fulfilment systems with low utilisation, provided capacity is not used efficiently and thus the order fulfilment system is not competitive in the long run due to the recently intensified competition (Van Gils et al. 2018; Kundu et al. 2020). On the other hand, in order fulfilment systems with a static bottleneck, the most crucial measure for improvement is not workload balancing but increasing the processing performance at the bottleneck. Consequently, we expect order fulfilment systems with high utilisation and a shifting bottleneck to be the main application field for workload balancing. For performance analysis of these systems, one should prefer the exact model as the simplified model is limited to worst-case analyses of system throughput and service level.

## **8.3 Memory and Computation Time Requirements**

In the simplified model, any multi-stage order fulfilment system is modelled as a single-stage system consisting of one aggregated process (see Figure 5.2). Thus, modelling the corresponding Markov chain is less complex than modelling the one in the exact model. The simplified model can be seen as a special case of the exact model since its Markov chain ensues from the Markov chain of the exact model when modelling is restricted to single-stage order fulfilment systems. Consequently, we expect differences between the simplified and the exact model regarding memory usage and computation time. In the following, we initially derive two hypotheses on the impact of the selected model for performance analysis on memory and computation time requirements (see Section 8.3.1). Subsequently, we verify these hypotheses in a numerical analysis (see Section 8.3.2).

#### **8.3.1 Hypotheses**

In the simplified model, any multi-stage order fulfilment system is modelled as a single-stage system by combining the processing steps p ∈ P to one aggregated process. The process index p is neglected in the definition of system state (see Section 7.2.1). System state is a (R + emax + 1)-dimensional vector in the simplified model, whereas system state in the exact model corresponds to a (pmax × (R + emax + 1))-dimensional matrix (see Section 6.2.1). The size of the state space in the exact model (see equation (6.45)) can be expressed as an integer multiple N ∈ N of the size of the state space in the simplified model (see equation (7.9)) with

$$N = \prod\_{p=2}^{p\_{max}} \left[ \prod\_{k=-R}^{-1} \left( (e\_{max} + |k|) \cdot h\_{(p-1),max} + 1 \right) \right. \tag{8.29}$$

$$\cdot \prod\_{k=0}^{e\_{max}} \left( (e\_{max} - k + 1) \cdot h\_{(p-1),max} + 1 \right) \Big| \,. \tag{8.20}$$

Thus, the state space of the Markov chain in the exact model is always larger than the one in the simplified model. As the size of the state space drives the memory requirements to store the state space, the transition matrix, and the limiting distribution of the Markov chain as well as the computational effort to calculate the Markov chain (see Section 6.4), we formulate the following hypotheses regarding memory and computation time requirements:

Hypothesis 5. Performance analysis of order fulfilment systems based on the exact model requires more memory than performance analysis based on the simplified model.

Hypothesis 6. Performance analysis of order fulfilment systems based on the exact model requires more computation time than performance analysis based on the simplified model.

### **8.3.2 Numerical Results**

To verify hypotheses 5 and 6, we conduct a numerical analysis based on samples C and E. Sample C contains 256 two-stage order fulfilment systems, and sample E consists of 353 three-stage order fulfilment systems. The samples differ regarding the ranges of the parameters, especially regarding the discretisation of the number of incoming orders per time period A: Its range is given by A<sup>C</sup> = {3, 4, . . . , 14} in sample C and A<sup>E</sup> = {2, 3, . . . , 7} in sample E. A comprehensive description of the samples is given in Appendix C.2. As evaluation criterion of the computation time, we measure the time required to calculate the Markov chain and the performance measures. Furthermore, we use the number of reachable states of the Markov chain as an indicator of memory requirements.

Table 8.6 presents the number of reachable states and the computation time of the exact and the simplified model as well as the absolute and relative deviations in the number of reachable states and the computation time between the exact and the simplified model for samples C and E. For each of these criteria, Table 8.6 gives the minimum, the maximum, and the average value.

In sample C, the number of reachable states is on average 135,072 states in the exact model and on average 7,568 states in the simplified model. In sample E, we observe on average 74,503 reachable states in the exact model and on average 1,101 reachable states in the simplified model. Thus, we observe an average relative deviation in the number of reachable states between the simplified and the exact model of -91.7% in sample C and -98.1% in sample E. These results confirm that performance analysis based on the exact model requires significantly more memory than performance analysis based on the simplified model (see hypothesis 5).

In sample C, the average computation time is 1.4 hours for the exact model and 2.6 seconds for the simplified model. In sample E, we observe an average computation time of 25.6 minutes for the exact model and 0.2 seconds for the simplified model. In both samples, the average relative deviation in computation time between the simplified and the exact model is -99%. These results confirm


Table 8.6: Number of reachable states and computation time of the exact and the simplified model.1,2

<sup>1</sup> We apply the strategies for runtime and memory optimisation introduced in Chapter 10 to calculate the Markov chain of the exact and the simplified model.

<sup>2</sup> Computations are conducted on a server with a CPU of 64 kernels and 128 threads and a RAM of 128 GB.

<sup>3</sup> Difference between the value of the considered criterion in the simplified model and the one in the exact model.

<sup>4</sup> Difference between the value of the considered criterion in the simplified model and the one in the exact model, divided by the value in the exact model.

that performance analysis based on the exact model requires significantly more computation time than performance analysis based on the simplified model (see hypothesis 6).

Note that it is not meaningful to compare the absolute values of the number of reachable states and the computation time between samples C and E since they depend on the discretisation of the stochastic parameters of the order fulfilment systems, which differs between samples C and E (see Appendix C.2).

## **8.4 Chapter Conclusion**

In this chapter, we evaluated and compared the exact and the simplified model for performance analysis regarding modelling accuracy, accuracy of performance analysis, memory usage, and computational effort.

Regarding memory usage and computational effort, the Markov chain of the simplified model requires significantly less memory and significantly shorter computation times than the corresponding Markov chain of the exact model due to its smaller number of reachable states (average reduction in the number of reachable states by more than 90%). Thus, when memory and computation time are scarce, the simplified model should be preferred over the exact model.

Regarding modelling accuracy, we showed that the simplified model suffers from modelling inaccuracies for performance analysis of order fulfilment systems with a shifting bottleneck and high utilisation: We mathematically proved that for this configuration of order fulfilment systems, the expected value of the system throughput calculated based on the simplified model is smaller than the actual expected value of the system throughput calculated based on the exact model. Furthermore, we numerically showed that these modelling inaccuracies of the simplified model lead to inaccuracies in the values of customer-related performance measures. For order fulfilment systems with high utilisation and a shifting bottleneck, the simplified model results in smaller values of α- and β-service level than the exact model. For instance, for systems with a utilisation of [0.9,1.0), the absolute value of the relative deviation of α-service level is at most 3.3% and on average 0.5% in the sample of two-stage systems, and it is at most 47.6% and on average 6.9% in the sample of three-stage systems. Hence, for order fulfilment systems with high utilisation and a shifting bottleneck, it is only reasonable to use the simplified model for worst-case analyses of system throughput and service level. In contrast, for order fulfilment systems with high utilisation and a static bottleneck, we mathematically proved that the expected value of the system throughput calculated based on the simplified model corresponds to the actual expected value of the system throughput calculated based on the exact model. For order fulfilment systems with low utilisation, we numerically showed that system performance calculated based on the simplified model is the same as the one calculated based on the exact model, apart from negligible numerical errors.

Order fulfilment systems with high utilisation and a shifting bottleneck are predominant in practical applications of workload balancing. For performance analysis of these systems, the exact model outperforms the simplified model in terms of modelling accuracy and accuracy of performance analysis, despite its drawbacks regarding memory and computation time requirements. The simplified model is limited to a worst-case performance analysis for these systems.

# **9 Formalisation and Solution Algorithms of the Capacity Planning Problem**

The analytical models developed in Chapters 6 and 7 enable the performance analysis of any order fulfilment system with a given system configuration. However, the focus of operations managers is not on the performance analysis of a given system configuration but on adapting the current system configuration to guarantee promised performance requirements to their customers. In particular, they have to decide on the capacity that is provided at the processing steps of the order fulfilment system to meet the performance requirements of the customers. Hence, this chapter aims at formalising the capacity planning problem in order fulfilment systems and at providing suitable solution algorithms. Section 9.1 introduces the mathematical formulation of the capacity planning problem, and Section 9.2 provides an overview of the research field of derivative-free and blackbox optimisation algorithms. In Section 9.3, we identify Mesh Adaptive Direct Search (MADS) and Surrogate Optimisation Integer (SO-I) to be suitable solution algorithms for capacity planning in order fulfilment systems. Their procedure and their problem-specific configuration are presented in Sections 9.4 and 9.5, respectively. By summarising the results of this chapter, Section 9.6 provides an answer to the third research question of this thesis.

## **9.1 Capacity Planning Problem**

Capacity planning in multi-stage order fulfilment systems determines the minimum possible capacity that has to be provided at each processing step to meet the performance requirements of the customers, such as a β-service level of 98%, in a multi-stage order fulfilment system with levelled order release. In the following, we introduce the mathematical formulation of the capacity planning problem (see Section 9.1.1) and specify its characteristics (see Section 9.1.2).

#### **9.1.1 Mathematical Formulation**

The decision variables of the capacity planning problem are defined by

$$\mathbf{c} = \begin{pmatrix} c\_1 & \dots & c\_{p\max} \end{pmatrix},\tag{9.1}$$

whereby c<sup>p</sup> ∈ N specifies the capacity provided at processing step p ∈ P. The objective function minimises the sum of provided capacity:

$$\min \sum\_{p \in \mathcal{P}} c\_p.$$

We ensure that the performance SL(c) achieved with the capacity c is as least as high as the required performance SL<sup>∗</sup> by formulating the following constraint:

$$SL(\mathbf{c}) \ge SL^\*.$$

Performance requirements in order fulfilment are commonly given as service level requirements. However, any other performance measure of the order fulfilment system (see Table 6.1) can be used to specify the performance requirements of the customers. It is even possible to consider multiple different performance requirements in the capacity planning problem. In this case, each performance requirement is modelled as a separate constraint.

Furthermore, we ensure the order fulfilment system to be stable by limiting the order income-related utilisation U˜(c) of the order fulfilment system with the provided capacity c to at most one:

$$
\vec{U}(\mathbf{c}) \le 1.
$$

We can derive lower cp,min and upper limits cp,max of the provided capacity cp, p ∈ P, based on the ranges of the number of incoming orders per time period A and the processing performances per time unit L<sup>p</sup> of the processes p ∈ P of the considered order fulfilment system. We obtain the minimum possible value of the required capacity c1,min at process p = 1 if the number of incoming orders per time period is at its minimum amin, and the processing performance per time unit of process p = 1 is at its maximum l1,max. For processes p ∈ P \ {1}, the minimum possible value of the required capacity cp,min corresponds to the smallest value of the minimum possible number of incoming orders per time period amin and the minimum possible processing performances per time period (cp<sup>0</sup> ,min · lp<sup>0</sup> ,min) of the previous processing steps (p <sup>0</sup> < p):

$$c\_{p,min} = \begin{cases} \left\lceil \frac{a\_{min}}{l\_{1,max}} \right\rceil & p = 1\\ \left\lceil \frac{\min\{a\_{min,\min\{p'\_{p'} \in \mathcal{P}, p'$$

In contrast, we obtain the maximum possible value of the required capacity c1,max at process p = 1 if the number of incoming orders per time period is at its maximum amax, and the processing performance per time unit of process p = 1 is at its minimum l1,min. For processes p ∈ P \ {1}, the maximum possible value of the required capacity cp,max results from the ratio of the maximum possible processing performance per time period c(p−1),max · l(p−1),max of the previous processing step (p − 1) to the minimum possible processing performance per time unit lp,min of process p:

$$c\_{p,max} = \begin{cases} \left\lceil \frac{a\_{max}}{l\_{1,min}} \right\rceil & p = 1\\ \left\lceil \frac{c\_{(p-1),max} \cdot l\_{(p-1),max}}{l\_{p,min}} \right\rceil & p \in \mathcal{P} \nmid \{1\}. \end{cases} \tag{9.3}$$

Consequently, the domain C of the decision variables c is specified as follows

$$\mathcal{C} = \{ \mathbf{c} \in \mathbb{N}^{p \times max} \mid c\_{p, min} \le c\_p \le c\_{p, max} \,\forall p \in \mathcal{P} \}\,. \tag{9.4}$$

The resulting mathematical formulation of the capacity planning problem is given as follows

$$\min \sum\_{p \in \mathcal{P}} c\_p \tag{9.5}$$

$$\text{s.t. } SL(\mathbf{c}) \ge SL^\* \tag{9.6}$$

$$
\bar{U}(\mathbf{c}) \le 1 \tag{9.7}
$$

$$\mathbf{c} \in \mathcal{C} \tag{9.8}$$

c ∈ N pmax . (9.9)

#### **9.1.2 Characteristics**

The capacity planning problem (9.5)-(9.9) is a constrained integer optimisation model. Its special characteristic is that the relationship between the provided capacity c and the performance SL(c) that is achieved with this capacity cannot be specified by a mathematical equation describing the performance of the order fulfilment system depending on the provided capacity. Instead, the relationship between the provided capacity and the resulting performance is specified by the analytical model for performance analysis that calculates for any order fulfilment system with given capacity the resulting system performance. This kind of optimisation model is called blackbox optimisation model.

The optimisation model

$$\min \; \Theta(\xi) \tag{9.10}$$

$$\text{s.t. } \Lambda\_j(\xi) \le 0 \qquad \forall j \in \mathcal{J} \tag{9.11}$$

	- ξ ∈ R n (9.13)

with objective function Θ(ξ), a finite number of constraints Λ<sup>j</sup> (ξ), j ∈ J , and a finite variable domain Ω is said to be a blackbox optimisation model if the objective function Θ(ξ) and/or the constraints Λ<sup>j</sup> (ξ), j ∈ J , are not given as mathematical equations but as a blackbox (Audet and Hare 2017, p.5).

In the case of the capacity planning problem, constraint (9.6) is given as a blackbox, whereby the analytical model for performance analysis represents the blackbox. In conclusion, the blackbox optimisation model for capacity planning has the following characteristics:


Furthermore, the blackbox constraint is said to be relaxable since it does not need to be satisfied to obtain a meaningful output of the blackbox (Le Digabel and Wild 2015). These characteristics become crucial in Section 9.3 when selecting suitable solution algorithms for the capacity planning problem.

## **9.2 Derivative-free and Blackbox Optimisation Algorithms**

In contrast to classical optimisation models, gradient information of the objective functions and the constraints of blackbox optimisation models are not or only partially available. Thus, common gradient-based solution approaches of optimisation models cannot be applied to solve blackbox optimisation models, such as the capacity planning problem. The research field of derivative-free and blackbox optimisation studies optimisation algorithms that do not require derivative information. Thus, these optimisation algorithms are applicable to solve the capacity planning problem.

## **9.2.1 Definition**

Derivative-free optimisation is defined as the mathematical study of optimisation algorithms that do not use gradient information. In contrast, blackbox optimisation focuses on the design and analysis of algorithms that assume the objective function and/or the constraints to be given as blackboxes. Blackbox optimisation usually does not make any assumptions regarding continuity, differentiability, or smoothness of the outputs. It includes research on heuristic methods or ad hoc methods for solving blackbox optimisation models. In contrast, derivative-free optimisation focuses on methods that can be mathematically analysed to prove convergence and/or to provide a stopping criterion (Audet and Hare 2017, p.5f.).

Typical application fields of derivative-free and blackbox optimisation are optimisation models whose functions are provided by a computer simulation that cannot be easily subjected to differentiation, optimisation models that involve conducting laboratory experiments that cannot be described by explicit mathematical equations, and optimisation models including noisy functions whose gradients are highly unreliable (Audet and Hare 2017, p.11f.).

If gradient information of the objective function and the constraints are available, reliable, and obtainable at reasonable cost, derivative-free and blackbox optimisation algorithms rarely outperform standard gradient-based optimisation algorithms. Thus, in this case, gradient-based algorithms should always be preferred (Audet and Hare 2017, p.6).

#### **9.2.2 Classification**

There are several criteria to classify the field of derivative-free and blackbox optimisation algorithms. Rios and Sahinidis (2013) provide a review on derivativefree algorithms for blackbox optimisation models with bound constraints and classify the algorithms into direct versus model-based algorithms, local versus global algorithms, and stochastic versus deterministic algorithms. The classification into direct and model-based methods is the most common one in the literature.

#### **9.2.2.1 Direct Search Methods**

There is no exact definition of the term direct search in the literature. It was first mentioned by Hooke and Jeeves (1961), who defined direct search as "sequential examination of trial solutions involving comparison of each trial solution with the best [solution] obtained up to that time together with a strategy for determining what the next trial solution will be." A direct search method is an iterative method that only uses comparisons of function values to determine trial points and that does not attempt to develop or use approximate gradients (Audet 2014; Larson et al. 2019, p.293). The pioneers among the direct search methods are the Coordinate Search (Fermi 1952), the Nelder-Mead algorithm (Nelder and Mead 1965), and the Hooke and Jeeves algorithm (Hooke and Jeeves 1961).

#### **9.2.2.2 Model-based Derivative-free Methods**

Model-based derivative-free methods iteratively construct a surrogate model to approximate the blackbox model (Sóbester et al. 2012). The objective function and the constraints of a surrogate model are reasonably accurate representations of the ones of the blackbox model at least locally, but the surrogate model is much faster to evaluate than the blackbox model. Note that the surrogate model of the blackbox optimisation model (9.10)-(9.13) consists of (|J |+ 1) surrogate functions, one for the objective function and one for every constraint j ∈ J . In each iteration, the surrogate model is used to determine new candidate points (Koziel et al. 2011, p.34f.), (Audet and Hare 2017, p.236).

The generic procedure of a model-based derivative-free method consists of the following phases (Vu et al. 2017):


In the design phase, a set of starting points is created, and the blackbox model is evaluated at these points. The Latin hypercube design is a commonly used approach to create this set. The purpose of the design phase is to create a set of points uniformly spread over the domain Ω of the blackbox model to get a first, global picture of its behaviour. In the model phase, a surrogate model of the blackbox model is constructed based on the already evaluated points. The most commonly used surrogate functions to approximate the blackbox model are polynomials, radial basis functions, and kriging models. It can be reasonable to consider multiple surrogate models in the model phase to prevent the user from choosing a poorly fitted model (Viana et al. 2013). Subsequently, in the search phase, a new candidate point is selected based on the surrogate model, the blackbox model is evaluated at this point, and this point is added to the set of already evaluated points. The selection of a new candidate point is the most crucial step in this phase. There is always the trade-off between exploitation – select a point around the current best point – and exploration – select a point that is far from all previously considered points to ensure that no global solution is overlooked. The model and the search phase are repeated until a predefined stopping criterion, such as a maximum number of evaluations of the blackbox model, is met (Vu et al. 2017).

Koziel et al. (2011) and Vu et al. (2017) provide comprehensive reviews of model-based derivative-free algorithms for blackbox optimisation models.

## **9.3 Selection of Solution Algorithms for Capacity Planning**

To select suitable algorithms for the capacity planning problem, we refer to the following criteria:


We select one direct search method and one model-based derivative-free method to investigate a representative of each category of blackbox optimisation algorithms regarding its performance in solving the capacity planning problem.

The focus of derivative-free and blackbox optimisation algorithms in the literature is on unconstrained problems on the one hand and optimisation models with continuous decision variables on the other hand. However, derivative-free algorithms solving constrained integer blackbox optimisation models have barely been studied in the literature so far (Müller et al. 2014).

Regarding the class of direct search methods, we prefer the Mesh Adaptive Direct Search (MADS) (Audet and Dennis Jr 2006) over the Generalised Pattern Search (Torczon 1997), and the Coordinate Search (Fermi 1952) since the MADS algorithm is the only method that considers general constraints. Furthermore, integer variables can be easily handled by the MADS algorithm since its mesh imposes a discrete structure on the variable domain. The MADS algorithm is a local optimisation method whose convergence to a local optimum is mathematically proven (Audet 2014). The implementation of the MADS algorithm is provided in the open-source software tool NOMAD (Abramson et al. 2018).

Regarding the class of model-based derivative-free methods, we select the Surrogate Optimisation Integer (SO-I) algorithm (Müller et al. 2014), which is an asymptotically complete model-based algorithm for integer optimisation models with computationally expensive blackbox objective functions and constraints. The implementation of the SO-I algorithm is available as MATLAB surrogate model toolbox MATSuMoTo (Müller 2014).

## **9.4 Capacity Planning using Mesh Adaptive Direct Search (MADS)**

The direct search method Mesh Adaptive Direct Search (MADS) introduced by Audet and Dennis Jr (2006) is an iterative procedure to identify a local optimum of the blackbox optimisation model by evaluating the blackbox model at selected points, starting at a given initial point. In each iteration, new points are selected based on different criteria: Search step and poll step. The main characteristic of the MADS algorithm is that every trial point has to be part of the mesh. The mesh is a discretisation of the variable domain Ω, and the mesh size parameter controls its fineness (Audet and Hare 2017, p.136-140).

In the following, we initially provide a high-level description of the MADS algorithm for solving unconstrained blackbox optimisation models with bound constraints (see Section 9.4.1). Subsequently, we present several selected extensions of the MADS algorithm, such as constraint handling, that are relevant for solving the capacity planning problem (see Section 9.4.2). Finally, we introduce the problem-specific configuration of the MADS algorithm for capacity planning in order fulfilment systems (see Section 9.4.3).

## **9.4.1 General Procedure**

We only provide a high-level description of the procedure and the main characteristics of the MADS algorithm for solving unconstrained blackbox optimisation models with bound constraints. A detailed description can be found in Audet and Hare (2017, p.136-140).

Each iteration of the MADS algorithm consists of the following steps:


In the search step, any method can be used to select a finite number of trial points. The focus of the search step is on exploring the variable domain Ω to escape from local optima and seek a global optimum. Line search, surrogate models, or heuristic methods, such as Tabu search, variable neighbourhood search, and simulated annealing, are commonly used methods (see Section 9.4.2.4). The set of trial points is evaluated either opportunistically or completely. In the case of an opportunistic evaluation, the evaluation of trial points stops as soon as a new incumbent solution is found. Otherwise, in the case of a complete evaluation, all trial points are evaluated, and the best trial point is selected as new incumbent solution, only if it outperforms the current incumbent solution.

In the poll step, a finite number of trial points around the current incumbent solution is selected based on given polling conditions:


The MADS algorithm provides different strategies to generate the polling directions (see Section 9.4.2.2). The set of trial points is evaluated using either the opportunistic or the complete strategy, analogous to the search step. The poll step ensures global convergence towards a local optimum.

In the step of parameter update, mesh and poll size parameter are updated depending on whether the current iteration was successful – a new incumbent solution was found – or unsuccessful – no new incumbent solution was found. The mesh size parameter converges much faster towards zero than the poll size parameter does. Hence, trial points in the search and the poll step can be chosen on a finer mesh than the mesh defined by the poll size parameter.

At the end of each iteration, the stopping criterion is checked. Either the total number of blackbox model evaluations, the computation time, or a threshold of the mesh size parameter are used to specify the stopping criterion. If the stopping criterion is met, the algorithm terminates. Otherwise, a new iteration is conducted.

The choice of the initial point is a crucial issue since the number of required iterations significantly reduces when choosing a good initial point. Either the initial point is chosen based on practitioners' professional intuition, or a small number of blackbox model evaluations is conducted to explore the behaviour of the blackbox model in the variable domain and to determine the initial point based on this sample.

### **9.4.2 Extensions**

In the following, we consider several extensions of the generic procedure of the MADS algorithm presented in Section 9.4.1. However, we limit ourselves to the extensions that are relevant for solving the capacity planning problem by the MADS algorithm.

#### **9.4.2.1 Constraint Handling**

The MADS algorithm provides multiple approaches for constraint handling depending on the fact whether the constraints are relaxable or unrelaxable. A relaxable constraint does not need to be satisfied to obtain a meaningful output of the blackbox. In contrast, a blackbox model with an unrelaxable constraint only provides a meaningful output if the unrelaxable constraint is satisfied (Le Digabel and Wild 2015).

The extreme barrier approach considers unrelaxable constraints. It sets the objective value to infinity for every infeasible point. Thus, the constrained blackbox optimisation model is treated as an unconstrained one, and any feasible point is preferred over any infeasible one (Audet and Dennis Jr 2006).

The progressive barrier approach uses the following non-negative constraint violation function

$$\Gamma(\xi) = \sum\_{j \in \mathcal{J}} (\max\left\{0; \Lambda\_j(\xi)\right\})^2 \tag{9.14}$$

to handle relaxable constraints. It places a threshold on the constraint violation it allows. All infeasible points that exceed this threshold are rejected, whereby the value of the threshold is progressively tightened with an increasing number of conducted iterations. The MADS algorithm with progressive barrier approach differentiates between feasible and infeasible incumbent solutions. In the poll step, polling is done around some points of the set of feasible incumbent solutions and the set of infeasible incumbent solutions (Audet and Dennis Jr 2009).

The progressive to extreme barrier approach is a hybrid method of progressive and extreme barrier approach. Relaxable constraints are initially treated by the progressive barrier approach. If polling around an infeasible incumbent solution generates a new infeasible incumbent that satisfies a constraint violated by the current infeasible incumbent, the corresponding constraint is treated by the extreme barrier approach in all future iterations. Consequently, all relaxable constraints are treated by the extreme barrier approach after finitely many iterations (Audet et al. 2010).

### **9.4.2.2 Strategies to Generate Poll Directions**

The MADS algorithm provides several strategies to generate the directions used to construct trial points in the poll step. The precondition is that the directions form a positive spanning set.

GPS is the most straightforward strategy to generate poll directions. Following the approach of the Generalised Pattern Search (Torczon 1997), GPS uses the coordinate directions as poll directions.

LTMADS is a random procedure to generate poll directions. The main drawbacks are that the runs of the MADS algorithm are not reproducible due to the randomly generated directions and that the generated directions are not necessarily orthogonal, which possibly leads to large angles between directions. Audet and Dennis Jr (2006) provide a detailed description of LTMADS.

In contrast, ORTHOMADS is a deterministic procedure to generate poll directions that provides an orthogonal positive spanning set of polling directions. Thus, the directions cover the surface of the unit sphere more densely and evenly than the ones of LTMADS. The procedure of ORTHOMADS is presented in detail in Abramson et al. (2009).

#### **9.4.2.3 Surrogate Management Framework**

If it is computationally expensive to evaluate the blackbox model at a certain point, it can be reasonable to use surrogate models. An overview of commonly used surrogate functions is given in Table 9.1. The surrogate management framework specifies how surrogate models can be integrated into direct search methods, especially into the MADS algorithm (Audet 2014).

On the one hand, surrogate models can be applied to conduct a model-based search in the search step. For this, a surrogate model is constructed based on the set of already evaluated points. The surrogate model is optimised to generate a Table 9.1: Commonly used surrogate functions in derivative-free and blackbox optimisation (Conn and Le Digabel 2013; Vu et al. 2017; Audet et al. 2018; Bhosekar and Ierapetritou 2018).


finite set of trial points, and the blackbox model is evaluated at these trial points using the opportunistic strategy (Audet and Hare 2017, p.238).

On the other hand, surrogate models can be used to order a given set of trial points either in the search or the poll step. For this, a surrogate model is constructed based on the set of already evaluated points. The surrogate model is evaluated at every trial point, and the trial points are sorted according to their value of the surrogate model, such that the most promising trial point is the first point of the ordered set. Finally, the blackbox model is evaluated at the ordered trial points using the opportunistic strategy (Audet and Hare 2017, p.238).

### **9.4.2.4 Methods to Use in the Search Step**

There is a plethora of suitable methods to create trial points in the search step. It is even possible to use several different methods in the search step. In the following, we focus on the methods applied to solve the capacity planning problem by the MADS algorithm.

As already described in Section 9.4.2.3, a model-based method can be used in the search step to create a set of trial points based on a surrogate model.

Speculative search represents a further suitable method for the search step. If the previous iteration of the MADS algorithm succeeded in finding an improved incumbent solution, the speculative search in the search step of the current iteration executes a simple line search along the previously successful polling direction (Audet and Dennis Jr 2006).

Moreover, it is possible to use a problem-specific adaption of the Nelder-Mead algorithm in the search step. Each iteration of the adapted Nelder-Mead algorithm starts with a set of affinely independent points defining a simplex. These simplex points are ordered based on the values of objective function and constraints from the best to the worst. In each iteration, the worst point of the simplex is replaced by a new point. Audet and Tribes (2018) give a detailed description of the adapted Nelder-Mead algorithm.

## **9.4.3 Problem-specific Configuration**

As already indicated by the discussed extensions in Section 9.4.2, the MADS algorithm provides a plethora of algorithmic parameters enabling its problemspecific configuration. The user manual of the corresponding software tool NOMAD (Audet et al. 2009) provides a comprehensive overview of all algorithmic parameters.

Regarding a problem-specific configuration of the MADS algorithm for the capacity planning problem, we determine the most important algorithmic parameters based on the results of numerical studies provided in the literature: First, we choose the progressive barrier approach (see Section 9.4.2.1) to handle the blackbox constraint of the capacity planning problem since it is a relaxable constraint (see Section 9.1.2). Furthermore, numerical studies show that the progressive barrier approach is the most reliable strategy for constraint handling (Audet and Dennis Jr 2009; Audet et al. 2010). Second, regarding the strategy to generate the poll directions (see Section 9.4.2.2), numerical studies show that ORTHOMADS outperforms GPS and that solutions of ORTHOMADS are not significantly bested by LTMADS in any scenario, especially in the case of computationally expensive blackbox models, when only one run is conducted (Abramson et al. 2009; Audet et al. 2010). Since the blackbox of the capacity planning problem is computationally expensive, we prefer ORTHOMADS over LTMADS.

We use the speculative search in the search step (see Section 9.4.2.4) and we will investigate the added value of an additional use of the model-based method and the Nelder-Mead algorithm in a numerical analysis on fine-tuning of the MADS algorithm for solving the capacity planning problem in Chapter 11. Furthermore, we will investigate different approaches to determine the initial point of the MADS algorithm in this numerical analysis. The stopping criterion of the MADS algorithm is defined by a threshold on the mesh size parameter. We choose a tolerance of one since this corresponds to the natural threshold of the mesh size parameter concerning integer blackbox optimisation models.

The blackbox optimisation model for capacity planning (9.5)-(9.9) has the following characteristics:


• For every point within the variable domain C, the value of the objective function and the value of U˜(c) are finite, and the value of SL(c) is in the range [0,1] due to the domain of the service level (see Section 6.3.7).

Based on these characteristics, applying the MADS algorithm with the progressive barrier approach, the constraint violation function (9.14), and OR-THOMADS to the capacity planning problem generates a convergent refining subsequence of incumbent solutions. The limit of the convergent refining subsequence is a Clarke KKT stationary point for the optimisation of the capacity planning problem. Audet et al. (2010) provide details on the convergence analysis of the MADS algorithm with progressive barrier approach and OR-THOMADS.

## **9.5 Capacity Planning using Surrogate Optimisation Integer (SO-I)**

The Surrogate Optimisation Integer (SO-I) algorithm introduced by Müller et al. (2014) is a model-based derivative-free method for solving constrained, integer blackbox optimisation models. In the following, we provide a highlevel description of the SO-I algorithm (see Section 9.5.1) and introduce the problem-specific configuration of the SO-I algorithm for capacity planning in order fulfilment systems (see Section 9.5.2).

### **9.5.1 General Procedure**

The SO-I algorithm consists of two separate optimisation phases. The first optimisation phase aims at finding a first feasible point if the initial set of points does not contain any feasible point. For this, a model-based derivative-free method (see below) is used to minimise the constraint violation function

$$\Gamma(\xi) = \sum\_{j \in \mathcal{J}} \max\{0; \Lambda\_j(\xi)\} \tag{9.15}$$

with respect to the bound and integrality constraints of the blackbox optimisation model. The second optimisation phase aims at identifying the minimum of the constrained blackbox optimisation model. For this, a model-based derivativefree method (see below) is used to minimise the penalty-augmented objective function

$$\Theta\_p(\xi) = \begin{cases} \Theta\_{\max} + \rho \cdot \sum\_{j \in \mathcal{J}} \left( \max\{ 0; \Lambda\_j(\xi) \} \right)^2 & \text{if } \xi \text{ is infeasible} \\ \Theta(\xi) & \text{if } \xi \text{ is feasible,} \end{cases} \tag{9.16}$$

with respect to the bound and integrality constraints of the blackbox optimisation model. Θmax denotes the current worst feasible objective function value, and ρ denotes the penalty factor.

The model-based derivative-free method used in both optimisation phases follows the generic procedure of model-based derivative-free methods as described in Section 9.2.2.2: In the design phase, a Latin hypercube design is used to create the initial set of points. A cubic radial basis function with l2-norm and a linear polynomial tail is used as surrogate function. In the search phase, a finite set of trial points is created by


Based on a score, the new candidate point is selected from the set of trial points. The score is calculated as the weighted sum of the predicted objective function value on the one hand and a distance criterion on the other hand. The predicted objective function value of a trial point corresponds to the value of the surrogate model at this point, scaled to the range [0,1]. The distance criterion calculates the distance of a trial point to the set of already evaluated points, whereby the result is also scaled to the range [0,1]. The iterative procedure terminates if either a local minimum is found, or the number of blackbox model evaluations equals the predefined maximum number of blackbox model evaluations (Müller et al. 2014).

The SO-I algorithm is asymptotically complete. Under the assumption of indefinitely long runtime and exact calculations, the SO-I algorithm will find the global minimum of the blackbox optimisation model with probability one (Müller 2014).

## **9.5.2 Problem-specific Configuration**

The SO-I algorithm, as described in Section 9.5.1, does not contain any degree of freedom enabling a problem-specific configuration of the algorithm with respect to the capacity planning problem. However, the SO-I implementation in the MATLAB-toolbox MATSuMoTo provides several algorithmic parameters to configure the SO-I algorithm. We will investigate some of them, such as the selected surrogate function and the selected sampling strategy in the search step, in a numerical analysis in Chapter 11 in order to fine-tune the SO-I algorithm for solving the capacity planning problem.

## **9.6 Chapter Conclusion**

In this chapter, we formulated the decision problem of capacity planning in multistage order fulfilment systems with performance requirements as a mathematical optimisation model and provided two solution algorithms.

The capacity planning problem is a blackbox optimisation problem since the relationship between the provided capacity and the performance that is achieved with this capacity cannot be specified by a mathematical equation, but it is given by the analytical model for performance analysis. The analytical model for performance analysis that calculates for any order fulfilment system with given capacity the resulting system performance represents the blackbox of the capacity planning problem. The research field of derivative-free and blackbox optimisation provides solution algorithms for blackbox optimisation models that do not require derivative information. We selected the direct search method Mesh Adaptive Direct Search (MADS) and the model-based derivative-free method Surrogate Optimisation Integer (SO-I) as suitable solution algorithms for capacity planning since they meet the characteristics of the capacity planning problem, their convergence is mathematically proven, and their implementation is open source.

The problem-specific configurations of the MADS algorithm and the SO-I algorithm enable a target-oriented determination of the minimum required, process-specific capacity to meet any performance requirement of the customers that is specified based on one or multiple performance measures of the order fulfilment system (see Table 6.1). Thus, the MADS algorithm and the SO-I algorithm provide an answer to the third research question of the thesis:

How can we determine the capacity required to meet specific performance requirements in multi-stage, stochastic order fulfilment systems with levelled order release and customerrequired order deadlines?

# **10 Runtime and Memory Optimisation of the Markov Chain**

When modelling system behaviour of a real-life order fulfilment system as a discrete-time Markov chain, memory and computation time requirements become challenging issues. This chapter aims at evaluating multiple approaches for runtime and memory optimisation of a Markov chain. In Section 10.1, we introduce multiple strategies to optimise memory and computation time requirements of a Markov chain. The optimisation potentials of these strategies regarding the models provided in this thesis are analysed in a numerical study in Section 10.2. Section 10.3 summarises the results of this chapter.

## **10.1 Strategies for Runtime and Memory Optimisation**

Memory and computation time requirements of a Markov chain are driven by the size of its state space and the discretisation of its stochastic parameters (see Section 6.4). Since the size of the state space of the Markov chain quickly increases when modelling the system behaviour of larger order fulfilment systems (see equation (6.45)), memory usage and computational effort are often challenging issues. In the following, we propose multiple strategies to optimise runtime and memory usage of a Markov chain.

## **10.1.1 Limitation of State Space**

The state space of the Markov chain results from its lower and upper limit (see Section 6.2.3), and its size is calculated by equation (6.45). However, depending on the parameterisation of the order fulfilment system and the discretisation of its stochastic parameters, this set of states contains multiple unreachable states. Unreachable states are not relevant in the subsequent calculation steps of the Markov chain. Consequently, these states are removed from the state space, and the subsequent calculation steps are limited to the set of reachable states.

The classical procedure to determine the set of reachable states of the Markov chain initially calculates all theoretically reachable states based on the lower and upper limit of the state space as well as the corresponding transition matrix. Unreachable states are then identified based on the transition matrix. State i ∈ X is said to be an unreachable state if the i th column of the transition matrix is a zero-vector. Finally, the identified unreachable states are removed from the state space and the transition matrix to obtain the set of reachable states and the corresponding transition matrix. The main drawback of this procedure is its inefficient use of runtime and memory: State space and transition matrix are first calculated and stored for all theoretically reachable states. Subsequently, unreachable states are removed.

The alternative procedure, which we propose to determine the set of reachable states of the Markov chain, iteratively identifies all states that are reachable from a given initial state of the Markov chain, either by a single-step transition or by an indirect transition via a finite number of other states. In this way, the alternative procedure determines the set of reachable states in a runtime- and memoryefficient manner since the state space and the transition matrix do not contain any unreachable state at any intermediate step of the calculation procedure.

For irreducible Markov chains, the alternative procedure does not provide any advantage compared to the classical procedure since the Markov chain consists of a unique communicating class. However, for reducible Markov chains, the alternative procedure does not determine the whole set of reachable states of the Markov chain but only the relevant subset of reachable states that contains all states that are reachable from the initial state. Thus, by using the alternative procedure, computation time and memory usage are reduced.

Based on the alternative procedure, it is impossible to calculate the matrix of limiting distributions of a reducible Markov chain. In contrast, we can only calculate the limiting distribution of the initial state. Hence, from a mathematical point of view, the alternative procedure is insufficient since it does not enable a complete description of the asymptotic behaviour of the Markov chain. However, from a practical perspective, it is sufficient to know the limiting distribution of the initial state of the Markov chain since the initial state describes the current state of the considered system.

### **10.1.2 Sparse Storage Schemes**

Transition matrices of Markov chains are usually large-sized but sparse matrices. A sparse matrix is a matrix with very few nonzero elements. The key idea of sparse storage schemes is to neglect the zero elements and to only store nonzero matrix entries. At the same time, common matrix operations still have to be possible (Saad 2003, p.73).

The compressed sparse row (CSR) format is the most popular sparse storage scheme. It uses three arrays Ev, Ec, and E<sup>r</sup> to store the sparse matrix E:


There are multiple variations of the CSR format, such as the compressed sparse column (CSC) format, which stores columns instead of rows (Saad 2003, p.89f.).

### **10.1.3 Indirect Solution Methods for Linear Systems**

The limiting distribution of the Markov chain is obtained by solving a set of linear equations based on the transition matrix of the Markov chain. Academic literature provides a plethora of methods to solve sets of linear equations. These methods are classified into direct and indirect solution methods. Direct methods modify the parameter matrix and use a fixed number of operations to exactly solve the set of linear equations. Common examples are the Gaussian elimination and Grassmann's algorithm (Bolch 2006, p.118f.). On the contrary, indirect methods iteratively calculate estimates starting with an initial estimate of the unknown solution. The sequence of estimates is expected to converge towards the solution eventually. The iterative procedure terminates if the estimates are sufficiently close to the exact solution. Hence, indirect methods only provide an approximation of the exact solution (Bolch 2006, p.132). Table 10.1 provides a comparison of direct and indirect solution methods for linear systems. We point out some differences in detail in the following.


Table 10.1: Comparison of direct and indirect solution methods for linear systems (Stewart 1994, p.61f.), (Bolch 2006, p.104,118f.,132).

Direct methods modify the entries of the parameter matrix, so that original zero entries can become nonzero entries. Consequently, round-off errors occur, and it is difficult to use sparse storage schemes. In contrast, indirect methods only use matrix multiplications that do not alter the parameter matrix: Original zero entries stay zero entries, and the sparsity of parameter matrix is preserved. Hence, there are no round-off errors, and sparse storage schemes are applicable. The convergence rate and the computation time of indirect methods depend on multiple factors, such as the quality of initial estimate, the required solution accuracy, the values of the parameter matrix, the chosen method, and the use of a preconditioning technique (Stewart 1994, p.61f.), (Bolch 2006, p.104,118f.,132). Preconditioning techniques are used to improve the efficiency and robustness of indirect methods. A preconditioner transforms the original set of linear equations into a set of linear equations that has the same solution but that is likely to be easier solved by an indirect method (Saad 2003, p.261f.). Saad (2003, p.283-345) provides a detailed overview of standard preconditioning techniques. In general, indirect methods are more efficient in space and time than direct methods. Hence, indirect methods are preferred over direct methods when solving large sets of linear systems (Bolch 2006, p.104).

In this thesis, we use the indirect method Generalized Minimal Residual Method (GMRES) that is a popular method for solving large linear systems. GMRES is an iterative Krylov subspace method to solve a nonsymmetric set of linear equations that minimises the norm of the residual vector over a Krylov subspace in each iteration (Saad and Schultz 1986). The procedure of GMRES is presented in detail in Stewart (1994, p.190-206) and Saad (2003, p.164-184).

### **10.1.4 Parallel Computing**

Parallel computing speeds up the calculation of a Markov chain since multiple calculation steps of the Markov chain can be conducted in parallel: First, when calculating the state space using the alternative procedure (see Section 10.1.1), we can determine the states reachable from a given state in a single-step transition in parallel for different states. Second, when calculating the transition matrix, it is possible to calculate each row in parallel, and for a given row, we can calculate all entries of this row in parallel. Finally, when calculating the limiting distribution, the matrix-vector operations, which are used to solve the set of linear equations, can be executed in parallel.

For the calculations of this thesis, we use the available processors of the CPU to conduct the aforementioned calculation steps in parallel. Graphics processing units (GPUs) provide further potentials to reduce computation time and memory usage. GPUs are widely used in high-performance computing, and they provide high computing power and large memory bandwidth. However, new programming challenges occur, such as adapting the parallel algorithms to the architecture of parallel computers with GPUs (Khodja et al. 2014). The calculation of Markov chains using GPUs is beyond the scope of this thesis.

## **10.2 Evaluation of Optimisation Potentials**

In the following, we evaluate the potentials of the proposed strategies for runtime and memory optimisation regarding the exact model for performance analysis in a numerical analysis based on samples B and D. Sample B contains 133 singlestage order fulfilment systems, and sample D consists of 256 two-stage order fulfilment systems. The samples differ regarding the ranges of the parameters, especially regarding the discretisation of the number of incoming orders per time period A: Its range is given by A<sup>B</sup> = {3, 4, . . . , 15} in sample B and A<sup>D</sup> = {2, 3, . . . , 6} in sample D. A comprehensive description of the samples is given in Appendix C.2.

### **10.2.1 Limitation of State Space**

The alternative procedure, proposed in Section 10.1.1 to calculate the set of reachable states of a Markov chain based on a given initial state, directly determines the set of reachable states. In contrast, the classical procedure initially calculates the set of all theoretically reachable states and subsequently reduces this set to the set of reachable states. To show and to quantify the benefit of the alternative procedure, we compare the number of theoretically reachable states, which is calculated based on equation (6.45), with the number of reachable states. Note that the initial state of each order fulfilment system is randomly selected from its set of theoretically reachable states.

Table 10.2 presents the number of theoretically reachable states and the number of reachable states, as well as the absolute and relative deviations between these two metrics for samples B and D. For each of these criteria, Table 10.2 gives the minimum, the maximum, and the average value.


Table 10.2: Number of theoretically reachable states, number of reachable states, and potentials of state space limitation when using the alternative procedure.

<sup>1</sup> Difference between the number of reachable states and the number of theoretically reachable states.

<sup>2</sup> Difference between the number of reachable states and the number of theoretically reachable states, divided by the number of theoretically reachable states.

For the single-stage order fulfilment systems of sample B, the number of theoretically reachable states is between 8,125 and 15,376 states, whereas the number of reachable states is between 156 and 9,893 states. For the two-stage order fulfilment systems of sample D, the number of theoretically reachable states is between 692,055 and 8,471,463 states. In contrast, the number of reachable states is between 1,401 and 9,625 states. Thus, we observe an average reduction in the number of calculated states of 43.6% in sample B and 99.8% in sample D. These results show that the number of states calculated by the alternative procedure is significantly smaller than the one calculated by the classical procedure. Consequently, we recommend using the alternative procedure instead of the classical procedure to calculate the set of reachable states of the Markov chain in order to reduce memory usage and computational effort.

By comparing the magnitude of the number of theoretically reachable states between samples B and D, we observe that the number of processes of the considered order fulfilment system has a significant impact on the number of theoretically reachable states. The average number of theoretically reachable states is 14,115 states for the single-stage systems of sample B, whereas it is 4,534,143 states for the two-stage systems of sample D. The number of theoretically reachable states results from the lower and upper limit of the state space (see equation (6.45)). Since the number of processes pmax determines the size of the exponent in equation (6.45), it has a major impact on the number of theoretically reachable states (see Section 6.4). In contrast, its impact on the number of reachable states is significantly smaller. Consequently, we expect that the potential of the alternative procedure in optimising memory and computation time requirements increases with an increasing number of processes of the considered order fulfilment system.

Note that it is not meaningful to compare the absolute values of the number of reachable states between samples B and D since these depend on the discretisation of the stochastic parameters of the order fulfilment systems, which differs between samples B and D (see Appendix C.2).

### **10.2.2 Indirect Solution Methods for Linear Systems**

In Section 10.1.3, we proposed using an indirect solution method for linear systems to reduce the computational effort of the Markov chain. To quantify the savings in computation time resulting from using an indirect solution method instead of a direct method for solving linear systems, we compare the computation time required to calculate the exact model for performance analysis when using the Gaussian elimination with the one when using GMRES.

Table 10.3 presents the minimum, the maximum, and the average value of the computation time in samples B and D for the following specifications of the Markov chain: Use of the Gaussian elimination and sequential computing, use of GMRES and sequential computing, and use of GMRES and parallel computing. Furthermore, Table 10.3 gives the minimum, the maximum, and the average value of the reduction in computation time when using GMRES instead of the Gaussian elimination and when computing the Markov chain in parallel instead of sequential computation.

In sample B, when using the Gaussian elimination, computation time is at least 0.04 seconds, at most 432 seconds, and on average 220 seconds. In contrast, when using GMRES, computation time is between 0.03 and 308 seconds, and it is 121 seconds on average. Thus, by using GMRES instead of the Gaussian elimination, we reduce the computation time of the exact model by on average 46%. Sample B contains one data point for which computation time increases by 4% when using GMRES instead of the Gaussian elimination. In sample D, when using the Gaussian elimination, computation time is between 2 and 428 seconds, and it is 125 seconds on average. In contrast, when using GMRES, computation time is at least 1.7 seconds, at most 286 seconds, and on average 64 seconds. Hence, we observe an average reduction in the computation time of the exact model by 49% when using GMRES instead of the Gaussian elimination. These results confirm that using the indirect method GMRES to solve linear systems provides significant savings in computation time compared to using the Gaussian elimination. Note that it is not meaningful to compare the absolute values of the computation time between samples B and D since the computation time


Table 10.3: Computation time and savings of computational effort when using GMRES and parallel computing.<sup>1</sup>

<sup>1</sup> Computations are conducted on a server with a CPU of 64 kernels and 128 threads and a RAM of 128 GB.

<sup>2</sup> Relative deviation of the computation time of the exact model that uses GMRES and that is computed in sequential manner from the computation time of the exact model that uses the Gaussian elimination and that is computed in sequential manner.

<sup>3</sup> Relative deviation of the computation time of the exact model that uses GMRES and that is computed in parallel from the computation time of the exact model that uses GMRES and that is computed in sequential manner.

depends on the ranges of the parameters of the order fulfilment systems, which differ between samples B and D (see Appendix C.2).

Indirect methods, such as GMRES, suffer from a lack of result accuracy, whereas direct methods are said to calculate exact results. However, direct methods suffer from round-off errors (see Table 10.1). For this reason, we evaluate the result accuracy by comparing the values of selected performance measures – α-service level, β-service level, expected number of unprocessed orders, and expected number of completed orders – obtained when using GMRES with the ones obtained when using the Gaussian elimination. Table 10.4 presents the absolute deviations of the selected performance measures when using GMRES instead of the Gaussian elimination in samples B and D. For each performance measure, Table 10.4 gives the minimum, the maximum, and the average value, as well as the 2.5%- and the 97.5%-quantile of the absolute deviations.


Table 10.4: Deviation of selected performance measures when using GMRES instead of the Gaussian elimination.

<sup>1</sup> Difference between the value of the considered performance measure resulting when using GMRES and the one resulting when using the Gaussian elimination.

The absolute deviation of α-service level is between -1.77E-04 and 1.2E-03 in sample B and between -2.43E-04 and 1.11E-04 in sample D. The corresponding 95%-confidence interval is given by [-7.26E-05, 4.06E-05] in sample B and [-1.04E-04, 3.57E-05] in sample D. The absolute deviations of β-service level, the expected number of unprocessed orders, and the expected number of completed orders are of comparable magnitude. Note that the calculated deviations only quantify the deviations obtained when using two different methods to solve linear systems that both suffer from inaccuracies. The actual values of the performance measures are still unknown. In conclusion, we state the observed deviations to be negligible for all considered performance measures, and we recommend using GMRES instead of the Gaussian elimination to speed up the calculation of the Markov chain.

### **10.2.3 Parallel Computing**

In Section 10.1.4, we pointed out some calculation steps of the Markov chain that can be conducted in parallel. To quantify the savings in computation time of parallel computing, we compare the computation time required to calculate the exact model for performance analysis in sequential manner with the one when calculating the exact model in parallel.

In sample B, computation time when using parallel computing is between 0.05 and 6.5 seconds, and it is 3.5 seconds on average. Compared to the computational effort of sequential computing, we observe an average reduction in computation time of 95%. For one data point in sample B, we observe a higher computation time for parallel computing (0.05 seconds) than for sequential computing (0.03 seconds). In sample D, computation time when using parallel computing is between 0.3 and 4.7 seconds, and it is 2.6 seconds on average. In comparison to the computational effort of sequential computing, we observe an average reduction in computation time of 95% (see Table 10.3). The observed savings in computation time when calculating the Markov chain in parallel confirm the potential of parallel computing for runtime optimisation of the Markov chain.

The outlier in sample B, for which the computation time of parallel computing is higher than the one of sequential computing, indicates that parallel computing is only beneficial for sufficiently large problem instances. For small problem instances, sequential computing can be quicker than parallel computing since the time required for starting and accumulating the threads in parallel computing exceeds the time saved by parallel computing. In conclusion, we recommend using parallel computing to speed up the calculation of the Markov chain.

## **10.3 Chapter Conclusion**

Runtime optimisation and efficient memory usage become challenging issues when modelling system behaviour of real-life order fulfilment systems as a discrete-time Markov chain since the size of the state space of the Markov chain quickly increases. For this reason, we proposed and analysed different strategies for runtime and memory optimisation of the Markov chain in this chapter.

By the proposed alternative procedure to calculate the set of reachable states of a Markov chain based on a given initial state, we achieve an average reduction in number of calculated states by 43.6% for single-stage order fulfilment systems and 99.8% for two-stage order fulfilment systems. A further strategy to reduce memory usage of the Markov chain is to use a sparse storage scheme to store the transition matrix of the Markov chain since transition matrices are usually sparse.

To reduce the computation time of the Markov chain, we recommend using the indirect method GMRES to solve linear systems and computing the Markov chain in parallel. By using GMRES instead of the Gaussian elimination, we achieve an average reduction in computation time by 46% for single-stage order fulfilment systems and 49% for two-stage order fulfilment systems. An analysis regarding the result accuracy when using GMRES indicates that the absolute deviations of the values of the performance measures are negligible. By calculating the Markov chain in parallel, we achieve an average reduction in computation time of 95% compared to a sequential computation. Parallel computing is limited to the available processors of the CPU. The additional use of GPUs is beyond the scope of this thesis. However, we expect the calculation of the Markov chain to further speed up when using GPUs.

# **11 Evaluation of Capacity Planning Algorithms**

This chapter aims at fine-tuning and evaluating the solution algorithms proposed in Chapter 9 for capacity planning in order fulfilment systems: Mesh Adaptive Direct Search (MADS) and Surrogate Optimisation Integer (SO-I). In Section 11.1, we specify the evaluation criteria. Both the MADS algorithm and the SO-I algorithm provide multiple algorithmic parameters for a problemspecific configuration of the algorithm. In Sections 11.2 and 11.3, we evaluate and compare several parameter settings of both algorithms to identify the most suitable parameter setting for the capacity planning problem, respectively. In Section 11.4, we evaluate and compare the overall performance of the algorithms. Based on the results of this chapter, Section 11.5 provides recommendations which algorithm to use for capacity planning in order fulfilment systems.

## **11.1 Evaluation Criteria**

We use solution quality and runtime efficiency as evaluation criteria for the numerical analysis and comparison of both different parameter settings of an algorithm and different algorithms based on a given sample of data points in the subsequent sections.

### **11.1.1 Solution Quality**

We evaluate the solution quality of an algorithm based on the objective function value of its solution (see equation (9.5)). Since the true minimum objective function value is unknown, we use the smallest objective function value of the solutions of all considered algorithms as an estimate of the minimum objective function value. To quantify the solution quality of an algorithm, we calculate the following key figures:


### **11.1.2 Runtime Efficiency**

We evaluate the runtime efficiency of an algorithm based on the number of calculated blackbox instances since in the case of capacity planning, the calculation of a blackbox instance, which corresponds to the calculation of the analytical model for performance analysis, is the most computationally expensive step of the algorithm. To quantify the runtime efficiency of an algorithm, we calculate the following key figures:


initially calculate the rank of each algorithm regarding the absolute number of calculated blackbox instances per data point. Subsequently, we obtain the average rank of runtime efficiency by calculating the average value of the rank values for all data points.

The average rank of runtime efficiency is a meaningful additional key figure since it abstracts from the differences in the absolute number of calculated blackbox instances between the data points of the considered sample. The smaller the average number of blackbox instances and the smaller the average rank of runtime efficiency of an algorithm, the higher the runtime efficiency of this algorithm.

## **11.2 Fine-Tuning of the MADS Algorithm**

The MADS algorithm provides a plethora of algorithmic parameters enabling a problem-specific configuration of the algorithm. We can determine some of them based on the results of numerical studies provided in the literature (see Section 9.4.3). However, further algorithmic parameters require a closer examination in a problem-specific numerical analysis. In the following, we initially determine the algorithmic parameters to investigate in the fine-tuning of the MADS algorithm (see Section 11.2.1). Subsequently, we evaluate and compare the resulting parameter settings to identify the most suitable parameter setting of the MADS algorithm for the capacity planning problem (see Section 11.2.2).

#### **11.2.1 Investigated Algorithmic Parameters**

The fine-tuning of the MADS algorithm focuses on different approaches to choose the initial point and different methods to use in the search step. Regarding the choice of the initial point, we either calculate the initial point based on equation (4.1) proposed in the step of system parametrisation of the levelling concept (see Section 4.3.1), or we use the best point of a Latin hypercube design whose sample size corresponds to 25% of the total number of points within the variable domain C of the capacity planning problem. Furthermore, we investigate using a surrogate model both as an additional method in the search step and to order the set of trial points in the search and the poll step (see Section 9.4.2.3). We consider five different surrogate functions: A quadratic model, a polynomial response surface, a radial basis function, a kernel smoothing model, and a Kriging model. Moreover, we investigate using the Nelder-Mead algorithm as an additional method in the search step.

Table 11.1 provides an overview of the values of the algorithmic parameters to investigate in the fine-tuning of the MADS algorithm. In total, we consider 24 different parameter settings of the MADS algorithm.


Table 11.1: Algorithmic parameters and their values investigated in the fine-tuning of the MADS algorithm.

### **11.2.2 Numerical Results**

We conduct a numerical analysis based on sample G to evaluate and compare the parameter settings of the MADS algorithm shown in Table 11.1 regarding solution quality and runtime efficiency. Sample G consists of 200 two-stage order fulfilment systems. A comprehensive description of sample G is given in Appendix C.3. The performance requirement of the capacity planning problem is specified by a predefined, individual value of α-service level. The results of the numerical analysis are given in Table 11.2.

For parameter settings SP/SMn/NMn and SP/SMy-QM/NMy, we obtain one unsolvable data point, respectively. Apart from this, the considered parameter settings do not differ regarding solution quality. Independent of the considered parameter setting, the MADS algorithm determines the optimal solution for all data points of sample G.

However, we observe differences regarding runtime efficiency between the considered parameter settings of the MADS algorithm. The average number of blackbox instances is between 20.19 and 24.03 for parameter settings using SP, whereas it is between 30.23 and 31.94 for parameter settings using LHD. The average rank of runtime efficiency is between 4.33 and 9.56 for parameter settings using SP and between 15.69 and 17.69 for parameter settings using LHD. Thus, using a Latin hypercube design to determine the initial point has a negative impact on the runtime efficiency. In contrast, regarding the use of Nelder-Mead algorithm and the use of a surrogate model, we do not observe any systematic impact on the runtime efficiency. We observe SP/SMy-QM/NMn (20.19), SP/SMn/NMn (20.61), and SP/SMy-QM/NMy (20.76) to be the best parameter settings regarding the average number of blackbox instances. SP/SMy-PRS/NMn (4.33), SP/SMy-KM/NMn (4.34), and SP/SMy-QM/NMn (4.62) are the parameter settings with the smallest average rank of runtime efficiency. Furthermore, LHD/SMy-KM/NMn (31.71; 17.12), LHD/SMy-RBF/NMn (31.7; 17.22), and LHD/SMy-PRS/NMn (31.94; 17.69) are the parameter settings with the highest average number of blackbox instances and the highest average rank of runtime efficiency, respectively.


Table 11.2: Solution quality and runtime efficiency of different parameter settings of the MADS algorithm based on sample G.

These results indicate that concerning the capacity planning problem, parameter settings using SP should be preferred over the ones using LHD. In conclusion, we recommend the parameter setting SP/SMy-QM/NMn of the MADS algorithm for solving the capacity planning problem in order fulfilment systems since the optimal solution is found for all data points of sample G when using this parameter setting. Furthermore, this parameter setting is the only parameter setting that belongs to the three best parameter settings regarding both evaluation criteria of runtime efficiency.

## **11.3 Fine-Tuning of the SO-I Algorithm**

The implementation of the SO-I algorithm in the MATLAB-toolbox MAT-SuMoTo provides some algorithmic parameters for a problem-specific configuration of the SO-I algorithm. Compared to the configuration options of the MADS algorithm, the ones of the SO-I algorithm are rather limited. In the following, we initially determine the algorithmic parameters to investigate in the fine-tuning of the SO-I algorithm (see Section 11.3.1). Subsequently, we evaluate and compare the resulting parameter settings to identify the most suitable parameter setting of the SO-I algorithm for the capacity planning problem (see Section 11.3.2).

### **11.3.1 Investigated Algorithmic Parameters**

The SO-I algorithm uses a Latin hypercube design to create the initial set of points in the design phase. We determine its sample size by 25% of the total number of points within the variable domain C of the capacity planning problem. A predefined maximum number of blackbox model evaluations is one of the stopping criteria of the SO-I algorithm. We determine the maximum number of blackbox evaluations by the total number of points within the variable domain C of the capacity planning problem.

The fine-tuning of the SO-I algorithm focuses on different sampling strategies and different surrogate functions. Regarding the sampling strategy that determines the finite number of trial points in the search phase, we either perturb the best feasible point found so far and uniformly select integer points in the variable domain, or we select the local minimum of the surrogate model as a new candidate point. Regarding the surrogate function, we use either a cubic radial basis function or a quadratic polynomial response surface.

Table 11.3 provides an overview of the values of the algorithmic parameters to investigate in the problem-specific fine-tuning of the SO-I algorithm. In total, we consider four different parameter settings of the SO-I algorithm.


Table 11.3: Algorithmic parameters and their values investigated in the fine-tuning of the SO-I algorithm.

### **11.3.2 Numerical Results**

We conduct a numerical analysis based on sample G (see Appendix C.3) to evaluate and compare the parameter settings of the SO-I algorithm shown in Table 11.3 regarding solution quality and runtime efficiency. The results of the numerical analysis are given in Table 11.4.

We obtain the optimal solution for all data points of sample G when using the sampling strategy CANDglob, independent of the selected surrogate function. In contrast, when using the sampling strategy SurfMin, we obtain a suboptimal solution for some data points of sample G. The share of suboptimal data points is 24% for parameter setting SurfMin/PRS and 1.5% for parameter setting


Table 11.4: Solution quality and runtime efficiency of different parameter settings of the SO-I algorithm based on sample G.

SurfMin/RBF. For both parameter settings, the corresponding average deviation from optimum equals one. These results indicate that using the sampling strategy SurfMin has a negative impact on the solution quality of the SO-I algorithm when solving the capacity planning problem.

Regarding runtime efficiency, we observe SurfMin/PRS (58.32) and CANDglob/RBF (58.36) to be the best parameter settings regarding the average number of blackbox instances. CANDglob/RBF (1.78) and SurfMin/RBF (1.85) are the parameter settings with the smallest average rank of runtime efficiency.

In conclusion, we recommend not to use parameter settings SurfMin/PRS and SurfMin/RBF of the SO-I algorithm due to the occurrence of suboptimal data points. Instead, we state parameter setting CANDglob/RBF to be the most suitable parameter setting of the SO-I algorithm for capacity planning since it outperforms parameter setting CANDglob/PRS in both evaluation criteria of runtime efficiency.

## **11.4 Comparison of the MADS Algorithm and the SO-I Algorithm**

In the following, we evaluate and compare the overall performance of the MADS algorithm and the SO-I algorithm to solve the capacity planning problem in order fulfilment systems in a numerical analysis based on sample G (see Appendix C.3). For this, we focus on the most suitable parameter setting of both algorithms for solving the capacity planning problem, respectively: Parameter setting SP/SMy-QM/NMn of the MADS algorithm and parameter setting CANDglob/RBF of the SO-I algorithm. We use the complete enumeration as a benchmark for the algorithms. The complete enumeration identifies the optimal solution of the capacity planning problem by investigating every point in the variable domain C.

We concentrate on a comparison regarding runtime efficiency since the results of Sections 11.2 and 11.3 show that the MADS algorithm and the SO-I algorithm do not differ regarding solution quality. The selected parameter settings of both algorithms identify the optimal solution for every data point of sample G.

Table 11.5 presents the minimum, the maximum, and the average value of the number of calculated blackbox instances per data point of the MADS algorithm, the SO-I algorithm, and the complete enumeration in sample G. Furthermore, Table 11.5 gives the minimum, the maximum, and the average value of the absolute and the relative reduction in the number of calculated blackbox instances when using the MADS algorithm compared to the complete enumeration, the SO-I algorithm compared to the complete enumeration, and the MADS algorithm compared to the SO-I algorithm.

The number of calculated blackbox instances per data point varies between 9 and 42 blackbox instances when using the MADS algorithm and between 49 and 72 blackbox instances when using the SO-I algorithm. The average number of calculated blackbox instances per data point is 20.19 for the MADS algorithm and 58.36 for the SO-I algorithm. The complete enumeration calculates 72 blackbox instances for every data point since the variable domain C of the


Table 11.5: Comparison of the MADS algorithm and the SO-I algorithm with complete enumeration regarding runtime efficiency based on sample G.

capacity planning problem is given by C = {1, 2, . . . , 6} × {1, 2, . . . , 12} (see equation (9.4)) for every data point of sample G. Thus, the total number of points within the variable domain is 72 points.

By using the MADS algorithm instead of complete enumeration, we observe an average reduction in the number of calculated blackbox instances per data point by 51.82 blackbox instances. The corresponding relative reduction is between 42% and 88%, and it is 72% on average. By using the SO-I algorithm, the number of calculated blackbox instances per data point reduces by on average 13.65 blackbox instances compared to complete enumeration. The corresponding relative reduction is on average 19% and at most 32%. For 12% of the data points of sample G, the number of calculated blackbox instances per data point of the SO-I algorithm corresponds to the total number of points within the variable domain. Thus, for these data points, the SO-I algorithm provides no benefit in runtime efficiency compared to complete enumeration. Furthermore, we observe an average reduction in the number of calculated blackbox instances per data point by 38.2 blackbox instances when using the MADS algorithm instead the SO-I algorithm. The corresponding relative reduction is between 26% and 88%.

In conclusion, these results indicate that the MADS algorithm has a remarkable higher runtime efficiency compared to the SO-I algorithm and complete enumeration. Consequently, we recommend using the MADS algorithm to solve the capacity planning problem in order fulfilment systems.

## **11.5 Chapter Conclusion**

In this chapter, we investigated multiple settings of algorithmic parameters of the MADS algorithm and the SO-I algorithm regarding solution quality and runtime efficiency in a numerical analysis to identify the most suitable parameter setting of each algorithm for capacity planning in order fulfilment systems. Furthermore, we evaluated and compared the overall performance of the MADS algorithm and the SO-I algorithm regarding solution quality and runtime efficiency.

We identify the parameter settings SP/SMy-QM/NMn and CANDglob/RBF to be the most suitable parameter settings of the MADS algorithm and the SO-I algorithm for solving the capacity planning problem, respectively. The MADS algorithm and the SO-I algorithm do not differ regarding solution quality since their most suitable parameter setting finds the optimal solution for all considered data points. However, the MADS algorithm has a remarkable higher runtime efficiency than the SO-I algorithm (average reduction in the number of calculated blackbox instances per data point by 65%) and complete enumeration (average reduction in the number of calculated blackbox instances per data point by 72%). Consequently, we recommend using the parameter setting SP/SMy-QM/NMn of the MADS algorithm for solving the capacity planning problem.

# **12 Evaluation of the Strategy of Levelled Order Release**

This chapter aims at evaluating the Strategy of Levelled Order Release in multistage, stochastic order fulfilment systems with customer-required order deadlines by comparing its performance with the one of an alternative strategy. A suitable alternative strategy is selected in Section 12.1. In Section 12.2, we derive several hypotheses on the expected behaviour of the strategies to compare. The models developed in this thesis enable an evaluation of the Strategy of Levelled Order Release regarding resulting system performance (see Section 12.3) as well as regarding required capacity (see Section 12.4). Section 12.5 summarises the results of this chapter.

## **12.1 Alternative Strategies**

The main characteristics of the Strategy of Levelled Order Release are (1) that a fixed capacity is reserved for order processing in each time period, and (2) that the orders are processed according to ascending due dates in each time period (see Chapter 4). Thus, the Strategy of Levelled Order Release comprises two different planning problems: (1) Capacity planning and (2) order dispatching. It is difficult to find suitable alternative strategies to evaluate the Strategy of Levelled Order Release since related research fields only focus on one of these planning problems.

Furthermore, to ensure comparability, the alternative strategies used to evaluate the Strategy of Levelled Order Release have to rely on the same amount of available information as the Strategy of Levelled Order Release:


In the following, we investigate the research fields of dispatching (see Section 12.1.1), scheduling (see Section 12.1.2), and flexible capacity adaption (see Section 12.1.3) to identify suitable alternative strategies for evaluating the Strategy of Levelled Order Release. The selected strategies are specified in Section 12.1.4.

### **12.1.1 Dispatching Policies**

A dispatching policy specifies a rule used to select the next order to be processed from a set of orders awaiting service (Blackstone et al. 1982). Blackstone et al. (1982), Pinedo (2009, p.442-444), and Baumann (2019, p.22-24) provide comprehensive reviews on dispatching policies in production systems.

Due to the limited focus of dispatching policies, the comparison of the Strategy of Levelled Order Release with any dispatching policy is limited to the planning problem of order dispatching. A comparison with respect to capacity planning is not possible. Instead, any dispatching policy selected to evaluate the Strategy of Levelled Order Release has to be combined with the capacity planning approach of the Strategy of Levelled Order Release in order to ensure comparability.

Dispatching policies are classified into the following categories:


Processing-related dispatching policies select the next order to be processed based on either the number of processing steps or the processing time. These policies are further differentiated regarding the considered metric – total sum or remaining sum – and the used sort order – ascending or descending. Common examples of this category are Smallest number of remaining operations, Highest number of total operations, Shortest remaining processing time, and Longest total processing time. The category of due-date related dispatching policies contains the policies Earliest due date and Minimum slack time. The slack time of an order specifies the time period that remains until reaching its due date and that is not required for its processing. The third category summarises further dispatching policies that use other characteristics to select the next order to be processed. Common examples of this category are Random service, First come first serve, and Greatest dollar first (Baumann 2019, p.22-24).

Processing-related dispatching policies implicitly assume each order to have a priori known, deterministic, and individual processing times and an individual processing sequence. In contrast, in order fulfilment systems, the processing performance per time period is stochastic, and it is specified by an order typeand processing step-specific probability distribution (see Section 3.2.3). Furthermore, all orders of a particular order type have the same processing sequence (see Section 3.2.2). Thus, orders of the same order type differ neither regarding processing time nor regarding processing sequence in order fulfilment systems. Since we investigate different order types separately in this thesis (see Section 5.3.1), both the processing time and the number of processing steps are no suitable criteria to select the next order to be processed in the order fulfilment system. Consequently, processing-related dispatching policies are unsuitable alternative strategies to evaluate the Strategy of Levelled Order Release.

The due date-related dispatching policy Earliest due date selects the next order to be processed based on its due date, analogous to the Strategy of Levelled Order Release. The dispatching policy Minimum slack time corresponds to Earliest due date in the context of order fulfilment due to the lack of order-specific processing times and an order-specific processing sequence. Consequently, both dispatching policies are unsuitable alternative strategies to evaluate the Strategy of Levelled Order Release.

Regarding the third category of dispatching policies, Greatest dollar first is an unsuitable alternative strategy since we do not model revenue and cost aspects of order fulfilment systems. However, Random service and First come first serve are suitable alternative strategies to evaluate the Strategy of Levelled Order Release.

## **12.1.2 Scheduling**

Scheduling is a decision-making process used in many manufacturing and service industries that allocates resources, such as machines and workers, to orders over a given planning horizon to optimise one or multiple objectives (Pinedo 2016, p.1). Scheduling is conducted based on the following information: A given finite set of orders, a given finite set of resources, the processing time of each order at each resource, and the release date, the due date, and the importance factor of each order (Pinedo 2016, p.13f.). In deterministic scheduling, processing times, release times, and due dates are deterministic, whereas, in stochastic scheduling, these parameters are stochastic, and scheduling is conducted based on the characteristics of their probability distributions (Pinedo 2016, p.246).

The comparison of the Strategy of Levelled Order Release with any scheduling problem is limited to the planning problem of order dispatching. A comparison with respect to capacity planning is not possible. Instead, any scheduling problem selected to evaluate the Strategy of Levelled Order Release has to be combined with the capacity planning approach of the Strategy of Levelled Order Release in order to ensure comparability.

In stochastic scheduling, processing time, release date, and due date of each order are defined by individual probability distributions, respectively. In contrast, in order fulfilment systems, the processing performance per time period is specified by an order type- and processing step-specific probability distribution (see Section 3.2.3), and the lead time is defined by an order type-specific probability distribution (see Section 3.2.2). Thus, orders of the same order type differ neither regarding processing times nor regarding lead time. Since we investigate different order types separately in this thesis (see Section 5.3.1), both processing time and lead time are unsuitable criteria to schedule the orders in an order fulfilment system. Consequently, stochastic scheduling is an unsuitable alternative strategy to evaluate the Strategy of Levelled Order Release.

### **12.1.3 Strategy of Flexible Capacity Adaption**

The Strategy of Levelled Order Release is characterised by a constant capacity per time period. In contrast, the main idea of the Strategy of Flexible Capacity Adaption is to match provided capacity as precisely as possible to demanded capacity resulting from the current system workload by hiring and laying-off workers (Chen et al. 2009).

There are different possible configurations of the Strategy of Flexible Capacity Adaption regarding the degree of capacity flexibility. In the case of an instantaneous capacity adaption, capacity levels are adapted at the beginning of each time period when order income and processing performance of the current time period are already known. Hence, provided capacity always fits exactly to demanded capacity. However, an instantaneous capacity adaption does not meet the assumptions of the Strategy of Levelled Order Release that order income and processing performance per time period are unknown in advance. Thus, it is an unsuitable alternative strategy to evaluate the Strategy of Levelled Order Release. Furthermore, periodic capacity adaptions are more realistic than instantaneous capacity adaptions for several reasons: First, working time decisions (working over/under time) are often taken periodically to abide to labour regulations and to accomplish the timely communication of these decisions to the relevant workers. Second, deployment of temporary workers may be restricted to specific times, such as the start of a day or the start of a week (Buyukkaramikli et al. 2013).

Independent of the chosen degree of flexibility, the Strategy of Flexible Capacity Adaption concentrates on capacity planning. The second planning problem of the Strategy of Levelled Order Release on order dispatching is not considered. Furthermore, the Strategy of Flexible Capacity Adaption leads to additional operational costs for hiring and laying off permanent and temporary workers and for other organisational efforts of capacity adaption. Cost aspects are neglected in our models of order fulfilment systems, but the evaluation of the Strategy of Levelled Order Release is restricted to the performance measures of an order fulfilment system (see Table 6.1). Consequently, in the context of this thesis, the Strategy of Flexible Capacity Adaption is an unsuitable alternative strategy since its flexibility costs cannot be incorporated into the comparison.

## **12.1.4 Selected Alternative Strategy**

Although none of the related research fields provides any planning approach that incorporates both capacity planning and order dispatching, we identify the dispatching strategies Random Service and First come first serve as suitable alternative strategies for evaluating the Strategy of Levelled Order Release. Random Service randomly selects the next order to be processed. First come first serve selects the next order to be processed based on the arrival times of the orders.

In the following, we focus on First come first serve (FCFS) to evaluate the Strategy of Levelled Order Release (LOR) since FCFS is widely used in production and logistics systems. Note that this comparison is limited to the planning problem of order dispatching. FCFS is combined with the capacity planning approach of LOR in order to ensure the comparability of both strategies.

## **12.2 Hypotheses**

The models developed in this thesis provide two possibilities to compare LOR and FCFS: When comparing the strategies regarding the resulting system performance, we investigate the impact of the selected strategy on the system performance, measured by at least one performance measure (see Table 6.1), for order fulfilment systems with a given capacity. Otherwise, when comparing the strategies regarding the required capacity, we analyse the impact of the selected strategy on the total required capacity in order fulfilment systems with predefined performance requirements.

In order fulfilment systems with low utilisation (see equation (8.1)), the order income per time period is the limiting factor of system throughput. The majority of incoming orders per time period is processed immediately after their time of arrival without being buffered at any processing step since there is enough idle capacity in systems with low utilisation. If orders are predominantly processed immediately after their time of arrival, their processing sequence has no impact on system throughput and the characteristics of the completed orders. These findings lead to the following hypothesis:

Hypothesis 7. In order fulfilment systems with low utilisation, system performance is independent of the strategy that is used to select the next order to be processed.

In contrast, in order fulfilment systems with high utilisation (see equation (8.2)), the processing performance per time period is the limiting factor of system throughput. A significant proportion of the incoming orders per time period is not processed immediately after their time of arrival, but the orders are buffered at some processing step of the multi-stage order fulfilment system before being processed. If it is impossible to process all incoming orders immediately after their time of arrival, a strategy is required to select the next order to be processed from the set of buffered orders. Thus, this strategy has an impact on system throughput and the characteristics of the completed orders.

FCFS does not consider the due dates of the orders, but the next order to be processed is selected based on its time of arrival. In contrast, LOR selects the order with the shortest due date as the next order to be processed. Hence, LOR tries to systematically avoid the occurrence of backorders. Its proportion of ontime processed orders is expected to be higher than the one of FCFS. α-service level measures the probability that none of the completed orders has a backlog duration, and β-service level quantifies the proportion of on-time completed orders on the total number of outgoing orders (see Section 6.3.7). Hence, we formulate the following hypothesis:

Hypothesis 8. In order fulfilment systems with high utilisation, the values of αand β-service level achieved when using LOR are at least as high as the ones achieved when using FCFS.

If it is possible to achieve higher values of α- and β-service level by using LOR instead of FCFS, as stated in hypothesis 8, this means the following in reverse with respect to capacity planning:

Hypothesis 9. The total required capacity to guarantee predefined performance requirements concerning α- and β-service level in order fulfilment systems when using LOR is at most as high as the one when using FCFS.

## **12.3 Comparison regarding Resulting System Performance**

In the following, we conduct a numerical analysis based on samples A1, A2, C, and E to compare LOR and FCFS regarding the resulting system performance and to verify hypotheses 7 and 8. We focus on α- and β-service level for simplicity reasons. As evaluation criteria, we calculate the relative deviation of the values of α- and β-service level achieved when using LOR from the ones achieved when using FCFS, respectively.

We use a simulation model to model the performance analysis in multi-stage, stochastic order fulfilment systems with FCFS. The simulation model is based on the one used to verify the exact model for performance analysis (see Appendix B). The simulation models only differ regarding the modelling of order processing in the simulation iteration when the next order to be processed is selected based on its time of arrival and not based on its due date. β-service level is used to determine the length of every simulation replication (precision of 1.0E-05) and the number of simulation replications (precision of 1.0E-04).

Samples A1 and A2 contain 120 and 350 single-stage order fulfilment systems, respectively. The ranges of the parameters are the same in both samples. Sample A2 focuses on order fulfilment systems with a utilisation of [0.5,1.0) since hypothesis 8 expects systems with high utilisation to be more relevant for the evaluation of LOR. Sample C contains 256 two-stage order fulfilment systems, and sample E consists of 353 three-stage order fulfilment systems. The samples differ regarding the ranges of the parameters, especially regarding the discretisation of the number of incoming orders per time period A: Its range is given by A<sup>A</sup> = {5, 6, . . . , 40} in samples A1 and A2, A<sup>C</sup> = {3, 4, . . . , 14} in sample C, and A<sup>E</sup> = {2, 3, . . . , 6} in sample E. A comprehensive description of the samples is given in Appendix C.2.

### **12.3.1 Analysis of** α**-Service Level**

Table 12.1 presents the relative deviations of α-service level between LOR and FCFS in samples A1, A2, C, and E. Each sample is subdivided into multiple classes according to the system utilisation (see equation (6.7)). Table 12.1 gives the minimum, the maximum, the average value, the 5%-, and the 95%-quantile of the relative deviations for each class of utilisation as well as for the whole sample.


Table 12.1: Deviation of α-service level between the Strategy of Levelled Order Release (LOR) and FCFS.


<sup>1</sup> Difference between the value of α-service level when using LOR and the one when using FCFS, divided by the value of α-service level when using FCFS.

We observe three data points each in sample A1 and E, four data points in sample C, and nine data points in sample A2 to have a negative relative deviation. The maximum of the negative relative deviations is 0.01% in sample A1, 0.02% in samples A2 and E, and 0.07% in sample C. The corresponding maxima of the negative absolute deviations are smaller than 5.3E-04. Thus, we state the observed negative deviations to be negligible. They result from numerical errors during the calculation. Thus, for all data points of the considered samples, α-service level of LOR is at least as high as the one of FCFS. Furthermore, we observe a maximum relative deviation of 141% in sample A1, 125% in sample A2, 74% in sample C, and 48% in sample E. Thus, there are data points in every sample for which α-service level achieved when using LOR is significantly higher than the one achieved when using FCFS.

Figure 12.1 presents the relative deviations of α-service level between LOR and FCFS depending on the utilisation of the order fulfilment system in samples A1, A2, C, and E. It indicates that the utilisation of the order fulfilment system has a

Figure 12.1: Deviation of α-service level between the Strategy of Levelled Order Release (LOR) and FCFS.

systematic impact on the relative deviations of α-service level since the relative deviations of α-service level increase as the utilisation increases.

In every sample, we observe a maximum relative deviation of at most 1.1% and an average relative deviation of at most 0.2% for systems with a utilisation smaller than 0.6. Consequently, for systems with a utilisation smaller than approximately 0.6, the deviations of α-service level between LOR and FCFS are negligible, and the strategy used to select the next order to be processed has a neglectable small impact on the resulting α-service level of the order fulfilment system.

In contrast, in sample A1, the relative deviation is 141% at its maximum and 54% on average for systems with a utilisation of [0.9,1.0). In sample A2, we observe a maximum relative deviation of 71% and an average relative deviation of 21% for systems with a utilisation of [0.9,0.95), whereas the relative deviation is 125% at its maximum and 51% on average for systems with a utilisation of [0.95,1.0). In sample C, we observe a maximum relative deviation of 74% and an average relative deviation of 21% for systems with a utilisation of [0.9,1.0). In sample E, the relative deviation is 48% at its maximum and 16% on average for systems with a utilisation of [0.9,1.0). These results indicate that for systems with high utilisation, the strategy used to select the next order to be processed has an impact on the resulting α-service level. For these systems, the α-service level achieved when using LOR can be significantly higher than the one when using FCFS.

By comparing the 90%-confidence intervals of the relative deviations of different classes of utilisation in each sample, we note that the width of the 90% confidence interval increases as the utilisation of the order fulfilment systems increases. In every sample, the 90%-confidence intervals per class of utilisation are narrower than 1% for systems with a utilisation smaller than 0.6. In contrast, for systems with a utilisation of [0.9,1.0), the 90%-confidence interval is given by, for instance, [1.8%, 137%] in sample A1 and [1%, 37%] in sample E. These results indicate that the variance of the relative deviations of α-service level increases as the utilisation of the order fulfilment system increases.

By comparing the magnitude of the relative deviations of α-service level between the samples, we observe that the number of processes of the considered order fulfilment system has a systematic impact on the relative deviation of α-service level. The maximum relative deviation is 141% for the single-stage systems of sample A1, 74% for the two-stage systems of sample C, and 48% for the three-stage systems of sample E. These results indicate that the benefit achieved when using LOR instead of FCFS decreases as the number of processes of the considered order fulfilment system increases.

### **12.3.2 Analysis of** β**-Service Level**

Table 12.2 presents the relative deviations of β-service level between LOR and FCFS in samples A1, A2, C, and E. Each sample is subdivided into multiple classes according to the system utilisation (see equation (6.7)). Table 12.2 gives the minimum, the maximum, the average value, the 5%-, and the 95%-quantile of the relative deviations for each class of utilisation as well as for the whole sample.

Relative deviation [%]<sup>1</sup> Utilisation Min Max Average Q0.05 Q0.95 No. Sample A1 [0.1,0.2) 0.00 0.00 0.00 0.00 0.00 4 [0.2,0.3) 0.00 0.00 0.00 0.00 0.00 14 [0.3,0.4) -0.01 0.04 0.00 -0.01 0.01 18 [0.4,0.5) 0.00 0.02 0.00 0.00 0.02 17 [0.5,0.6) 0.00 0.10 0.01 0.00 0.06 9 [0.6,0.7) 0.00 1.37 0.18 0.00 0.71 14 [0.7,0.8) -0.01 1.65 0.33 0.00 1.57 17 [0.8,0.9) 0.00 4.94 0.81 0.00 4.16 15 [0.9,1.0) 0.04 6.17 2.50 0.05 6.11 12 [0.1,1.0) -0.01 6.17 0.42 0.00 3.25 120 Sample A2 [0.5,0.55) -0.01 0.15 0.01 0.00 0.07 28 [0.55,0.6) -0.01 0.19 0.02 -0.01 0.14 39 [0.6,0.65) 0.00 0.48 0.05 0.00 0.28 40 [0.65,0.7) -0.01 1.17 0.12 0.00 0.65 34 [0.7,0.75) -0.01 1.16 0.19 0.00 0.94 37 [0.75,0.8) -0.01 1.99 0.33 0.00 1.80 49 [0.8,0.85) 0.00 4.49 1.13 0.00 3.31 35 [0.85,0.9) 0.00 5.71 1.06 0.00 4.02 32 [0.9,0.95) -1.78 6.26 1.29 -0.35 5.29 31

Table 12.2: Deviation of β-service level between the Strategy of Levelled Order Release (LOR) and FCFS.


<sup>1</sup> Difference between the value of β-service level when using LOR and the one when using FCFS, divided by the value of β-service level when using FCFS.

We observe a minimum relative deviation of -0.01% in all samples, except six data points in sample A2, whose relative deviation is between -7% and -0.7%, one data point in sample C with a relative deviation of -1.2%, and two data points in sample E with a relative deviation of -3.7%, respectively. Apart from these outliers that will be discussed in Section 12.3.3, we state the observed negative relative deviations to be negligible. They result from numerical errors during the calculation. Hence, we conclude that β-service level achieved when using LOR is at least as high as the one achieved when using FCFS. Furthermore, we observe a maximum relative deviation of approximately 6% in every sample. Thus, there are data points in every sample for which β-service level achieved when using LOR is significantly higher than the one achieved when using FCFS.

Figure 12.2 presents the relative deviations of β-service level between LOR and FCFS depending on the utilisation of the order fulfilment system in samples A1, A2, C, and E. It indicates that the utilisation of the order fulfilment system has a systematic impact on the relative deviations of β-service level since the relative deviations of β-service level increase as the utilisation increases.

Figure 12.2: Deviation of β-service level between the Strategy of Levelled Order Release (LOR) and FCFS.

In every sample, we observe a maximum relative deviation of at most 0.2% and an average relative deviation of at most 0.02% for systems with a utilisation smaller than 0.6. Consequently, for systems with a utilisation smaller than approximately 0.6, the deviations of β-service level between LOR and FCFS are negligible, and the strategy used to select the next order to be processed has a neglectable small impact on the resulting β-service level of the order fulfilment system.

In contrast, for systems with a utilisation of [0.9,1.0), the maximum relative deviation is between 6.1% and 6.7% in every sample. The corresponding average relative deviation is between 1.2% and 2.5%. These results indicate that for systems with high utilisation, the strategy used to select the next order to be processed has an impact on the resulting β-service level. In these systems, the β-service level achieved when using LOR can be significantly higher than the one achieved when using FCFS.

By comparing the 90%-confidence intervals of the relative deviations of different classes of utilisation in each sample, we note that the width of the 90% confidence interval increases as the utilisation of the order fulfilment systems increases. In every sample, the 90%-confidence intervals per class of utilisation are narrower than 1% for systems with a utilisation smaller than 0.7. In contrast, for systems with a utilisation of [0.9,1.0), the 90%-confidence interval is given by, for instance, [0.05%, 6.1%] in sample A1 and [0.07%, 5.4%] in sample E. These results indicate that the variance of the relative deviations of β-service level increases as the utilisation of the order fulfilment system increases.

### **12.3.3 Discussion**

For bothα- and β-service level, we do not observe significant differences between LOR and FCFS for order fulfilment systems with a utilisation smaller than approximately 0.6 in all samples. These results confirm that for order fulfilment systems with low utilisation, the strategy used to select the next order to be processed has a neglectable small impact on the resulting system performance (see hypothesis 7).

In contrast, in every sample, we observe significantly higher values of both αand β-service level when using LOR instead of FCFS in order fulfilment systems with high utilisation. For instance, for systems with a utilisation of [0.9,1.0), the average relative deviation over all samples is 27% for α-service level and 2% for β-service level. α-service level is more sensitive to the selected strategy than βservice level. These results confirm that the strategy used to select the next order to be processed has an essential impact on the resulting system performance in order fulfilment systems with high utilisation and that system performance achieved when using LOR can be significantly higher than the one achieved when using FCFS (see hypothesis 8).

Samples A2, C, and E contain few data points whose β-service level achieved when using LOR is smaller than the one achieved when using FCFS. Table 12.3 gives a detailed analysis of selected performance measures – order incomerelated utilisation, β-service level, and expected total number of lost sales per time period – of the corresponding order fulfilment systems. All of them are characterised by a high utilisation and a high proportion of backorders. In this kind of order fulfilment systems, LOR, which selects the next order to be processed based on the shortest due date to systematically avoid lost sales, has the negative side-effect that only a small proportion of the available processing performance per time period remains to process orders without failed due dates. Thus, the β-service level achieved when using LOR is relatively small in these systems. In contrast, FCFS, which does not rely on due dates when selecting the next order to be processed, accepts a higher number of lost sales, and thus more orders can be completed on time. Consequently, the β-service level achieved when using FCFS is higher than the one achieved when using LOR.


Table 12.3: Order income-related utilisation, β-service level, and expected total number of lost sales per time period of selected order fulfilment systems when using the Strategy of Levelled Order Release (LOR) and when using FCFS (A2.1-A2.6 of sample A2, C.1 of sample C, E.1-E.2 of sample E).

<sup>1</sup> Difference between the value of the key figure (SL<sup>β</sup> or E(S)) when using LOR and the one when using FCFS, divided by the value of the key figure when using FCFS.

## **12.4 Comparison regarding Required Capacity**

In the following, we conduct a numerical analysis based on samples F and G to compare LOR and FCFS regarding the total required capacity and to verify hypothesis 9. As evaluation criteria, we use the absolute and relative deviation of the total required capacity in the order fulfilment system to guarantee a predefined performance requirement when using LOR from the one when using FCFS. The total required capacity corresponds to the objective function value of the capacity planning problem (see equation (9.5)).

We use the parameter setting SP/SMy-QM/NMn of the MADS algorithm for solving the capacity planning problem (see Chapter 11). In the case of capacity planning for order fulfilment systems using LOR, the exact model for performance analysis represents the blackbox. In contrast, in the case of capacity planning for order fulfilment systems using FCFS, the adapted simulation model, which is used in Section 12.3 for performance analysis of order fulfilment systems using FCFS, represents the blackbox.

Sample F consists of 200 single-stage order fulfilment systems, and sample G contains 200 two-stage order fulfilment systems. The samples differ regarding the ranges of the parameters, especially regarding the discretisation of the number of incoming orders per time period A: Its range is given by A<sup>F</sup> = {3, 4, . . . , 15} in sample F and A<sup>G</sup> = {2, 3, . . . , 6} in sample G. A comprehensive description of the samples is given in Appendix C.3. In both samples, the performance requirement of the capacity planning problem is specified by a predefined, individual value of α-service level.

### **12.4.1 Numerical Results**

Table 12.4 presents the minimum, the maximum, and the average values of the absolute and relative deviations of the total required capacity between LOR and FCFS in samples F and G.


Table 12.4: Deviation of total required capacity between the Strategy of Levelled Order Release (LOR) and FCFS.

<sup>1</sup> Difference between the capacity required when using LOR and the one required when using FCFS.

<sup>2</sup> Difference between the capacity required when using LOR and the one required when using FCFS, divided by the capacity required when using FCFS.

We observe a maximum absolute deviation of zero in both samples. Thus, for all data points of the considered samples, the total required capacity to guarantee a predefined α-service level when using LOR is at most as high as the one when using FCFS. Furthermore, for 13.5% of the data points in sample F and 12% of the data points in sample G, the total required capacity to guarantee a predefined α-service level when using LOR is smaller than the one when using FCFS. We observe a minimum absolute deviation of -1.00 in each sample. The corresponding minimum relative deviation is -33% in sample F and -20% in sample G, respectively.

Focusing on the data points for which we observe a capacity saving when using LOR instead of FCFS, the average absolute capacity saving is one time unit in each sample. The corresponding average relative capacity saving is 21% in sample F and 14% in sample G (see Table 12.5).

Table 12.5: Capacity savings when using the Strategy of Levelled Order Release (LOR) instead of FCFS.


Figure 12.3 presents the relative deviation of the total required capacity between LOR and FCFS depending on the required α-service level in samples F and G. It indicates that the value of the required α-service level does not have any systematic impact on the capacity savings that can be achieved when using LOR instead of FCFS. The required α-service level of data points for which we observe a capacity saving when using LOR instead of FCFS is between 0.45 and 1.00 in sample F and between 0.08 and 1.00 in sample G.

### **12.4.2 Discussion**

We observe no data point with a positive deviation of total required capacity between LOR and FCFS in both samples. Furthermore, we observe capacity savings of on average 18% when using LOR instead of FCFS for approximately 13% of the data points in each sample. These results confirm that the total

Figure 12.3: Deviation of total required capacity between the Strategy of Levelled Order Release (LOR) and FCFS.

required capacity to guarantee a predefined performance requirement when using LOR is at most as high as the one when using FCFS (see hypothesis 9).

The absolute capacity saving equals one time unit for every data point, for which we observe a capacity saving when using LOR instead of FCFS. This is due to the fact that the considered order fulfilment systems in both samples are small (see Appendix C.3), and the capacity is an integer parameter (see equation (9.9)). Thus, the variable domain C of the capacity planning problem is limited: It is given by C <sup>F</sup> = {1, 2, . . . , 8} in sample F and by C <sup>G</sup> = {1, 2, . . . , 6} × {1, 2, . . . , 12} in sample G (see equation (9.4)).

## **12.5 Chapter Conclusion**

In this chapter, we evaluated the Strategy of Levelled Order Release in multistage, stochastic order fulfilment systems with customer-required order deadlines by comparing its performance with the one of FCFS.

The Strategy of Levelled Order Release comprises two planning problems: Capacity planning and order dispatching. It is difficult to find suitable alternative strategies to evaluate the Strategy of Levelled Order Release since related research fields only focus on one of these planning problems: Strategies of flexible capacity adaption concentrate on capacity planning, whereas dispatching strategies and scheduling problems focus on order dispatching. Consequently, the comparison of the Strategy of Levelled Order Release with any alternative strategy is always limited to one of its planning problems. The Strategy of Flexible Capacity Adaption is an unsuitable alternative strategy since the models provided in this thesis are restricted to performance analysis of order fulfilment systems and do not consider further evaluation criteria, such as the flexibility costs of the Strategy of Flexible Capacity Adaption. Furthermore, processing-related dispatching policies, such as Shortest remaining processing time, and scheduling problems are unsuitable alternative strategies due to modelling reasons, and duedate related dispatching policies select the next order to be processed based on the same criterion as the Strategy of Levelled Order Release. However, the dispatching strategy FCFS is a suitable alternative strategy to evaluate the Strategy of Levelled Order Release. Since FCFS is limited to the planning problem of order dispatching, it is combined with the capacity planning approach of the Strategy of Levelled Order Release in order to ensure the comparability of both strategies.

The comparison of the Strategy of Levelled Order Release and FCFS regarding the resulting system performance shows that the strategy used to select the next order to be processed has a neglectable small impact on the resulting values of α- and β-service level in order fulfilment systems with a utilisation smaller than approximately 0.6. In contrast, in order fulfilment systems with a utilisation higher than 0.6, one can achieve significantly higher values of α- and β-service level when using the Strategy of Levelled Order Release instead of FCFS. For instance, for systems with a utilisation of [0.9,1.0), we observe an average increase of α-service level of 27% and of β-service level of 2%.

The comparison of the Strategy of Levelled Order Release and FCFS regarding the total required capacity to guarantee a predefined α-service level in the order fulfilment system shows that the total required capacity when using the Strategy of Levelled Order Release is at most as high as the one when using FCFS. In approximately 13% of the considered capacity planning problems, we achieve a capacity saving of on average 18% when using the Strategy of Levelled Order Release instead of FCFS.

# **13 Conclusion**

This chapter aims at summarising the most important results of this thesis and presenting an outlook on further research topics.

## **13.1 Summary**

Order fulfilment systems are forced to efficiently manage a highly volatile customer demand consisting of low-volume orders while simultaneously meeting customer-required, short order deadlines. Hopp and Spearman (2004) provide three buffer types – inventory buffer, time buffer, and capacity buffer – to handle the volatile workload in production systems. In this thesis, we combine the potentials of time buffer and capacity buffer to develop and analytically investigate the Strategy of Levelled Order Release to balance workload in multi-stage, stochastic order fulfilment systems with customer-required order deadlines over time on a tactical level. The focus is on the two main application fields of order fulfilment: Warehouses and make-to-order systems. The contribution of this thesis is three-fold: We develop (1) a workload balancing concept, (2) a discrete-time analytical model for performance analysis, and (3) an algorithm for capacity planning under performance constraints for multi-stage, stochastic order fulfilment systems with levelled order release and customer-required order deadlines.

The results of the literature review state the research gap of this thesis. First, there is no workload balancing approach in the literature to balance the workload of order fulfilment systems over time on a tactical level by using a combination of time buffer and capacity buffer. Second, the discrete-time analytical model introduced in this thesis is the first analytical model for performance analysis of workload balancing in order fulfilment systems that ensures a real-life representation of order fulfilment systems. It regards customer-required order deadlines and the interdependencies between the processing steps, and it models customer demand, customer-required deadlines, and processing performance as stochastic parameters. Furthermore, the discrete-time analytical model enables the exact calculation of complete probability distributions for several systemand customer-related performance measures. Third, capacity planning problems with due date-related performance constraints that rely on probability distributions of corresponding performance measures and service level in multi-stage systems with stochastic, customer-required order deadlines have not been studied in the literature so far.

Workload Balancing in Order Fulfilment Systems The order fulfilment system forms the basis for all concepts and models developed in this thesis. We formally describe an order fulfilment system by a finite set of order types, a finite set of processes, and a function mapping each order type to its processing sequence in the order fulfilment system. Each order type is characterised by the number of incoming orders per time period and the lead time of an order, and each process is specified by its processing performance and its capacity. We develop a levelling concept for order fulfilment systems, the so-called Strategy of Levelled Order Release, based on the key ideas of Heijunka levelling in production systems. It is impossible to directly apply the concept of Heijunka levelling in order fulfilment systems due to several differences in the system characteristics, such as order lot size, the role of customer demand, and characteristics of capacity. Moreover, the overarching question of the levelling concept is different: Heijunka levelling decides on a suitable buffer size of the finishedgoods-inventory to meet the promised performance requirements, whereas, in order fulfilment systems, the decision is on a suitable amount of provided capacity to meet the promised performance requirements. The key characteristics of the Strategy of Levelled Order Release are the following:


In this way, the Strategy of Levelled Order Release systematically exploits the time buffer of each order between its time of arrival and its deadline to balance the variability of the customer demand. The remaining variability of the customer demand is either passed on to the customer, resulting in low service levels, or it is balanced using the capacity buffer. The extent to which capacity buffer is used to balance the remaining variability depends on the specific performance requirements of the order fulfilment system. Thus, the Strategy of Levelled Order Release provides an answer to the first research question of the thesis:

### How can we balance workload over time in multi-stage, stochastic order fulfilment systems with customer-required order deadlines?

Performance Analysis of Order Fulfilment Systems We choose a discretetime Markov chain to model and analyse system behaviour and performance of a multi-stage, stochastic order fulfilment system with levelled order release and customer-required order deadlines since discrete-time analytical models provide several advantages regarding accuracy, level of detail, and the description of real processes compared to continuous-time analytical models. Above all, a discrete-time Markov chain enables the exact calculation of complete probability distributions of stochastic performance measures. To model system behaviour of multi-stage order fulfilment systems with levelled order release as a discrete-time Markov chain, we consider the order types of the order fulfilment system separately. We assume each order type to have a sequential processing sequence, and we limit the possible backlog duration of an order by the maximum accepted backlog duration. Furthermore, regarding the modelling of partially

processed orders that are transmitted between the processing steps of a multistage order fulfilment system, we differentiate between an exact and a simplified modelling approach: The exact modelling approach provides an exact recording of partially processed orders, whereas partially processed orders are neglected in the simplified modelling approach. Based on these modelling approaches, we introduce an exact and a simplified analytical model for performance analysis of multi-stage, stochastic order fulfilment systems with levelled order release and customer-required order deadlines by


The simplified model can be seen as a special case of the exact model since it ensues from the exact model when modelling is restricted to single-stage order fulfilment systems. Based on the calculated performance measures, such as system throughput, α-, β- and γ-service level, utilisation, performance balance, number of unprocessed orders, number of lost sales, and time buffer and backlog duration of a completed order, the developed models enable a comprehensive performance analysis of multi-stage, stochastic order fulfilment systems with levelled order release and customer-required order deadlines. These models for performance analysis provide an answer to the second research question of the thesis:

### How can we determine the performance of multi-stage, stochastic order fulfilment systems with levelled order release and customerrequired order deadlines?

We show that the simplified model suffers from modelling inaccuracies for performance analysis of order fulfilment systems with high utilisation and a shifting bottleneck: We prove that for this system configuration, the expected value of system throughput calculated based on the simplified model is smaller than the actual one calculated based on the exact model. Furthermore, we numerically show that these modelling inaccuracies of the simplified model lead to inaccuracies in the values of α- and β-service level. For instance, for systems with a utilisation of [0.9,1.0), the absolute value of the relative deviation of α-service level is at most 3.3% and on average 0.5% in the sample of twostage systems, and it is at most 47.6% and on average 6.9% in the sample of three-stage systems. Hence, for order fulfilment systems with high utilisation and a shifting bottleneck, it is only reasonable to use the simplified model for worst-case analyses of system throughput and service level. In contrast, for order fulfilment systems with high utilisation and a static bottleneck, we prove that the expected value of system throughput calculated based on the simplified model corresponds to the actual one calculated based on the exact model. For order fulfilment systems with low utilisation, we numerically show that system performance calculated based on the simplified model is the same as the one calculated based on the exact model, apart from negligible numerical errors. However, we expect order fulfilment systems with high utilisation and a shifting bottleneck to be the main application field for workload balancing since systems with low utilisation are not competitive in the long run, and the most crucial measure for improvement in systems with a static bottleneck is not workload balancing but increasing the processing performance at the bottleneck. For performance analysis of these systems, one should prefer the exact model over the simplified model for reasons of modelling accuracy and accuracy of performance analysis, despite its drawbacks regarding computational effort and memory usage due to the larger number of reachable states (average reduction in the number of reachable states by more than 90% when using the simplified model instead of the exact one).

Capacity Planning in Order Fulfilment Systems The decision problem of capacity planning in multi-stage order fulfilment systems with performance requirements is a blackbox optimisation problem since the relationship between the provided capacity and the performance that is achieved with this capacity cannot be specified by a mathematical equation, but it is given by the analytical model for performance analysis. The analytical model that calculates the resulting system performance for any order fulfilment system with a given capacity represents the blackbox of the capacity planning problem. We select the direct search method Mesh Adaptive Direct Search (MADS) and the modelbased derivative-free method Surrogate Optimisation Integer (SO-I) as suitable solution algorithms for capacity planning since they meet the characteristics of the capacity planning problem, their convergence is mathematically proven, and their implementation is open source. The problem-specific configurations of the MADS algorithm and the SO-I algorithm enable a target-oriented determination of the minimum required, process-specific capacity to meet any performance requirement of the customers that is specified based on one or multiple performance measures of the order fulfilment system. Thus, both algorithms provide an answer to the third research question of the thesis:

How can we determine the capacity required to meet specific performance requirements in multi-stage, stochastic order fulfilment systems with levelled order release and customerrequired order deadlines?

Numerical Analysis and Application Runtime optimisation and efficient memory usage become challenging issues when modelling system behaviour of real-life order fulfilment systems as a discrete-time Markov chain since the size of the state space of the Markov chain quickly increases. Therefore, we introduce the alternative procedure to determine the set of reachable states of a Markov chain based on a given initial state. By this procedure, we achieve an average reduction in the number of calculated states by 43.6% for singlestage systems and by 99.8% for two-stage systems. A further strategy to reduce memory usage of the Markov chain is to use a sparse storage scheme to store the transition matrix of the Markov chain. To reduce the computation time of the Markov chain, we recommend using the indirect method GMRES to solve linear systems and computing the Markov chain in parallel. By using GMRES instead of the Gaussian elimination, we observe an average reduction in computation time by approximately 46%. By calculating the Markov chain in parallel, we achieve an average reduction in computation time by 95% compared to sequential computation.

A numerical analysis on fine-tuning and evaluating the solution algorithms for capacity planning regarding solution quality and runtime efficiency identifies the parameter settings SP/SMy-QM/NMn and CANDglob/RBF to be the most suitable parameter settings of the MADS algorithm and the SO-I algorithm for solving the capacity planning problem. The MADS algorithm and the SO-I algorithm do not differ regarding solution quality since their most suitable parameter settings find the optimal solution for all considered data points, respectively. However, the MADS algorithm has a remarkable higher runtime efficiency than the SO-I algorithm (average reduction in the number of calculated blackbox instances per data point by 65%). Consequently, we recommend using the parameter setting SP/SMy-QM/NMn of the MADS algorithm for solving the capacity planning problem in order fulfilment systems.

Finally, we evaluate the Strategy of Levelled Order Release in multi-stage, stochastic order fulfilment systems with customer-required order deadlines by comparing its performance with the one of FCFS. Since the Strategy of Levelled Order Release comprises two planning problems – capacity planning and order dispatching –, it is difficult to find suitable alternative strategies in related research fields to evaluate the Strategy of Levelled Order Release. Strategies of flexible capacity adaption concentrate on capacity planning, whereas dispatching strategies and scheduling problems focus on order dispatching. The Strategy of Flexible Capacity Adaption is an unsuitable alternative strategy since the models provided in this thesis are restricted to performance analysis of order fulfilment systems and do not consider further evaluation criteria, such as the flexibility costs of the Strategy of Flexible Capacity Adaption. Furthermore, processing-related dispatching policies, such as Shortest remaining processing time, and scheduling problems are unsuitable alternative strategies due to modelling reasons, and due-date related dispatching policies select the next order to be processed based on the same criterion as the Strategy of Levelled Order Release. However, the dispatching strategy FCFS is a suitable alternative strategy to evaluate the Strategy of Levelled Order Release. Since FCFS is limited to the planning problem of order dispatching, it is combined with the capacity planning approach of the Strategy of Levelled Order Release in order to ensure the comparability of both strategies. The comparison of the Strategy of Levelled Order Release and FCFS regarding the resulting system performance shows that the strategy used to select the next order to be processed has a neglectable small impact on the resulting values of α- and β-service level in order fulfilment systems with a utilisation smaller than approximately 0.6. In contrast, in order fulfilment systems with a utilisation higher than 0.6, one can achieve significantly higher values of α- and β-service level when using the Strategy of Levelled Order Release instead of FCFS. For instance, for systems with a utilisation of [0.9,1.0), we observe an average increase of α-service level of 27% and of β-service level of 2%. The comparison of the Strategy of Levelled Order Release and FCFS regarding the total required capacity to guarantee a predefined α-service level in the order fulfilment system shows that the total required capacity when using the Strategy of Levelled Order Release is at most as high as the one when using FCFS. In approximately 13% of the considered capacity planning problems, we achieve a capacity saving of on average 18% when using the Strategy of Levelled Order Release instead of FCFS.

## **13.2 Outlook**

After answering the research questions and summarising the results, we encounter the boundaries of this thesis. Some aspects go beyond the scope of this thesis but are interesting and relevant to analyse.

First, further research can be deducted on generalising the models developed in this thesis. In real-life order fulfilment systems, customer demand is often seasonally fluctuating. By modelling the number of incoming orders as a timedependent, stochastic parameter, the time-dependency of customer demand can be incorporated into the models of this thesis. Furthermore, the processing sequence of orders in order fulfilment systems usually contains multiple splits and merges, whereas the models of this thesis are limited to a sequential processing sequence. We can model an order type whose processing sequence contains splits by dividing this order type into multiple artificial order types, each of which has a sequential processing sequence. However, modelling merges in order fulfilment systems is beyond the scope of this thesis.

Second, parallel computing of the Markov chain is limited to the available processors of the CPU in this thesis. However, we expect the additional use of graphics processing units to provide further potentials in reducing computation time and memory usage of the Markov chain.

Finally, the concept of variability buffers of Hopp and Spearman (2004) provides further strategies for workload balancing in multi-stage, stochastic order fulfilment systems apart from the Strategy of Levelled Order Release. For instance, the Strategy of Flexible Capacity Adaption uses the capacity buffer to match provided capacity as precisely as possible to demanded capacity. Even inventory buffer provides some potentials for workload balancing in order fulfilment systems by doing preparatory processing steps in advance. Moreover, various combined strategies are reasonable. A comprehensive evaluation and comparison of these different approaches for workload balancing in multi-stage, stochastic order fulfilment systems require a holistic evaluation framework. Based on the models developed in this thesis, the comparison of the Strategy of Levelled Order Release with alternative strategies is limited to several performance measures so far. However, a holistic evaluation framework has to incorporate further evaluation criteria, such as cost aspects.

# **A Methodology of Literature Review**

The state of the art of academic literature in the research segments of this thesis, which is presented in Chapter 2, results from a keyword-based systematic literature review in the Scopus database. We searched for each pair of keywords given in Table A.1 in the article title, abstract, and keywords. We limited the resulting set of publications to journal articles and conference papers in English language. We additionally removed all papers of unrelated research fields. In the remaining relevant publications, we analysed the references and the citations to find further relevant publications that did not emerge previously.


Table A.1: Keywords of literature review.



# **B Simulation Model**

The simulation model depicts and analyses system behaviour and performance of multi-stage, stochastic order fulfilment systems with levelled order release and customer-required order deadlines. It is based on the exact modelling approach (see Figure 5.1), and it considers a multi-stage, stochastic order fulfilment system with a finite set of processes P, a unique order type, and levelled order release considering a finite set of due dates K. The order type and each process are specified by the parameters introduced in Sections 3.2.2 and 3.2.3, respectively. Analogous to the analytical models for performance analysis, we observe the order fulfilment system at discrete-time points in time t ∈ N0. The state of the simulation model specifies the number of unprocessed orders according to their due dates (see equation (6.8)).

## **B.1 Simulation Iteration**

A simulation iteration models order processing in the multi-stage order fulfilment system in any time period t. It consists of the same (pmax + 2) sub-steps as the state transition of the Markov chain of the exact model (see Section 6.2.2):

	- Order processing at process p = 1,
	- . . .,
	- Order processing at process pmax,

We use the same notation as in the exact model for performance analysis and denote the interim states of a simulation iteration by the variables y (m) , m ∈ {0, . . . ,(pmax + 2)}, whereby m is used as a counter, and y 0 specifies the state of the simulation before the start of order processing in time period t. The realisations of the order income per time period G = g and the processing performance per time period H = h in time period t are generated based on independent random generators. Depending on these realisations, the interim states y (m) , m ∈ {0, . . . ,(pmax + 2)}, of the simulation iteration in time period t are calculated using the equations of the state transition of the Markov chain of the exact model (see equations (6.9)-(6.12)). Based on the values of the interim states, we calculate the values of the performance measures of the order fulfilment system (see Table 6.1) in time period t as follows

• Number of unprocessed orders at process p at the beginning of the time period

$$q\_p = \sum\_{k \in \mathcal{K}} y\_{p,k}^{(0)} \qquad \forall p \in \mathcal{P} \tag{B.l}$$

• Number of unprocessed orders in the system at the beginning of the time period

$$q = \sum\_{p \in \mathcal{P}} \sum\_{k \in \mathcal{K}} y\_{p,k}^{(0)} \tag{B.2}$$

• Number of unprocessed orders at process p immediately before the start of order processing at process p

$$\bar{q}\_p = \sum\_{k \in \mathcal{K}} y\_{p,k}^{(p-1)} \qquad \forall p \in \mathcal{P} \tag{B.3}$$

• Number of lost sales per time period at process p

$$s\_p = y\_{p, -R}^{(p)} \qquad \forall p \in \mathcal{P} \tag{B.4}$$

• Total number of lost sales per time period

$$s = \sum\_{p \in \mathcal{P}} y\_{p, -R}^{(p)} \tag{B.5}$$

• Performance balance of process p

$$w\_p = h\_p - \sum\_{k \in \mathcal{K}} y\_{p,k}^{(p-1)} \qquad \forall p \in \mathcal{P} \tag{B.6}$$

• Order backlog-related utilisation of process p

$$u\_p = \min\left\{1; \frac{y\_{p,k}^{(p-1)}}{h\_p}\right\} \qquad \forall p \in \mathcal{P} \tag{B.7}$$

• Processed order per time period

$$\mathbf{m} = \begin{pmatrix} m\_{1, -R} & \dots & m\_{1, e\_{max}} \\ \vdots & \ddots & \vdots \\ m\_{p\_{max}, -R} & \dots & m\_{p\_{max}, e\_{max}} \end{pmatrix} \\ \tag{B.8}$$
 
$$m\_{p, k} = \min \left\{ y\_{p, k}^{(p-1)}; \max \left\{ 0; h\_p - \sum\_{l=-R}^{k-1} y\_{p, l}^{(p-1)} \right\} \right\} \\ \tag{B.8}$$
 
$$\forall p \in \mathcal{P}, k \in \mathcal{K}$$

• Number of processed orders per time period at process p

$$f\_p = \sum\_{k \in \mathcal{K}} m\_{p,k} \qquad \forall p \in \mathcal{P} \tag{B.9}$$

• Number of completed orders per time period

$$f = \sum\_{k \in \mathcal{K}} m\_{p\_{max},k} \tag{B.10}$$

• α-service level

$$sl\_{\alpha} = \begin{cases} 1 & \sum\_{k=-R}^{-1} m\_{p\_{max},k} = 0 \\ 0 & \text{otherwise} \end{cases} \tag{B.11}$$

• β-service level

$$sl\_{\beta} = \frac{\sum\_{k=0}^{e\_{max}} m\_{p\_{max},k}}{f+s} \tag{B.12}$$

• γ-service level

$$sl\_{\gamma} = 1 - \left(\frac{\sum\_{k=-R}^{-1} |k| \cdot m\_{pmax,k} + (R+1)\cdot s}{R \cdot f + (R+1)\cdot s}\right). \tag{B.13}$$

We register the time difference to the order deadline d, the backlog duration d backlog, and the time buffer d buf f er of all orders completed in time period t by recording the due date of every completed order at its time of completion.

### **B.2 Warm-up Period**

The simulation model suffers from the problem of initial transient. Since the initial state of the simulation is randomly selected, there is no guarantee that initial simulation iterations correctly represent system behaviour in steady-state. Instead, the average value of any key figure calculated based on a finite number of replications will be a biased estimator of the real value of this key figure in steady-state. To avoid this bias in estimators, a certain number of simulation iterations at the beginning of every replication, the so-called warm-up period of the simulation, is deleted, and the key figures are estimated based on the remaining observations that are said to represent steady-state behaviour of the simulation (Law 2015, p.511).

There are various methods in the literature to determine the end of the warm-up period τ . We use the Marginal Standard Error Rule (MSER-5) introduced by White Jr (1997) to identify the end of warm-up period τ based on β-service level and 100 replications, each of which consists of 5,000 simulation iterations.

## **B.3 Stopping Criteria**

The length of every replication and the number of replications are determined based on β-service level. SLβ,i,j denotes β-service level of simulation iteration j in replication i. SLβ,i denotes β-service level of replication i, and it is calculated as the average value of the values SLβ,i,j of β-service level of the simulation iterations in steady-state of replication i

$$SL\_{\beta,i} = \frac{1}{m\_i - \tau} \sum\_{k=\tau}^{m\_i} SL\_{\beta,i,k},\tag{B.14}$$

whereby m<sup>i</sup> denotes the number of simulation iterations of replication i.

To determine the length of replication i ∗ , we calculate the absolute deviation of the average value of β-service level measured after the current simulation iteration j ∗ from the average value of β-service level measured after simulation iteration j 0

$$\Delta\_{i^\*,j^\*,j'} = \left| \frac{1}{j^\* - \tau} \sum\_{k=\tau}^{j^\*} SL\_{\beta,i^\*,k} - \frac{1}{j' - \tau} \sum\_{k=\tau}^{j'} SL\_{\beta,i^\*,k} \right| \tag{B.15}$$

$$\forall j^\*, j' \ge \tau, (j^\* - j') \in \{1, 2, \dots, 10\}.$$

Replication i ∗ stops after simulation iteration j ∗ if the following condition holds

$$
\Delta\_{i^\*,j^\*,j'} < \epsilon\_R \qquad \forall (j^\*-j') \in \{1, 2, \dots, 10\}, \tag{B.16}
$$

whereby <sup>R</sup> denotes the predefined precision of a replication.

The replications of the simulation are independent and identically distributed since the initial state of every replication is selected randomly (Law 2015, p.489). To determine the number of replications, we calculate the 95%-confidence interval of β-service level based on the values SLβ,i of β-service level of the already calculated replications. The 95%-confidence interval of β-service level after replication i ∗ is calculated as follows

$$\begin{aligned} \left[ \overline{SL}\_{\beta,i^\*} - t\_{\{i^\*-1\},0.975}, \sqrt{\frac{S^2(SL\_{\beta,i^\*})}{i^\*}} \right] \\ \left[ \overline{SL}\_{\beta,i^\*} + t\_{\{i^\*-1\},0.975}, \sqrt{\frac{S^2(SL\_{\beta,i^\*})}{i^\*}} \right], \end{aligned} \tag{B.17}$$

whereby SLβ,i<sup>∗</sup> denotes the average value of β-service level after replication i ∗

$$\overline{SL}\_{\beta,i^\*} = \frac{1}{i^\*} \sum\_{k=1}^{i^\*} SL\_{\beta,k},$$

S 2 (SLβ,i<sup>∗</sup> ) denotes the variance of β-service level after replication i ∗

$$S^2(SL\_{\beta,i^\*}) = \frac{1}{i^\*-1} \sum\_{k=1}^{i^\*} (SL\_{\beta,k} - \overline{SL}\_{\beta,i^\*})^2,$$

and t(<sup>i</sup> <sup>∗</sup>−1),0.<sup>975</sup> denotes the 97.5%-quantile of the t-distribution with (i <sup>∗</sup> − 1) degrees of freedom. The simulation model terminates after replication i ∗ if the 95%-confidence interval of β-service level is narrower than the predefined precision <sup>S</sup> of the simulation

$$2 \cdot t\_{\{i^\*-1\}, 0.975} \cdot \sqrt{\frac{S^2(SL\_{\beta, i^\*})}{i^\*}} < \epsilon\_S. \tag{B.18}$$

## **B.4 Performance Measures of the Order Fulfilment System**

The values of each performance measure of the order fulfilment system (see Table 6.1) are derived from the corresponding absolute frequency distribution that results from the values of this performance measure in all steady-state simulation iterations of the simulation model. The probability distribution of each stochastic performance measure is estimated by the corresponding relative frequency distribution, whereas the value of each deterministic performance measure is estimated by the average value of the corresponding absolute frequency distribution.

# **C Design of Experiments**

We use a space-filling design to generate samples of order fulfilment systems. The main characteristic of space-filling designs is that the design becomes increasingly dense in the design space by increasing the sample size (Santner et al. 2018, p.149). Vazquez and Bect (2011) provide some justification for choosing space-filling designs by showing that a space-filling design is the design with the highest rate of the mean square prediction error decreases as design size increases. The Latin hypercube design, which we use in this thesis, is a popular space-filling design whose sample is evenly spread over the range of each input variable separately. Santner et al. (2018, p.152f.) provide a detailed description of the general procedure for obtaining a Latin hypercube sample. To determine a reasonable sample size, Jones et al. (1998) propose the rule of thumb to use a sample size of 10d data points for an input space of dimension d. Due to curse of dimensionality, a 10d sample becomes very sparse as d increases. Loeppky et al. (2009) investigate this issue and conclude that a sample size of 10d data points is reasonable for experiments with an input space of dimension d ≤ 5.

## **C.1 Latin Hypercube Designs for Order Fulfilment Systems**

An order fulfilment system with a unique order type is specified by the following parameters (see Sections 3.2 and 5.3):

• Number of incoming orders per time period A,


The maximum backlog duration R and the number of processes pmax are not included into the Latin hypercube design. They are both constant for every data point of a sample. The maximum backlog duration is given R = 1 in every sample. Regarding the number of processes pmax, we differentiate between samples of single-stage, two-stage, and three-stage order fulfilment systems. The remaining parameters are systematically varied in every sample.

Some of them – A, E, and Lp, p ∈ P – are stochastic parameters. To reduce the degree of freedom of the stochastic parameters and thus to reduce the number of input variables of the Latin hypercube design, we assume each stochastic parameter to belong to a particular class of theoretical probability distribution: We assume that the number of incoming orders per time period A is approximately Poisson-distributed with parameter λ and a finite range since arrival processes are commonly modelled as Poisson-distributed random variables in the literature (Law 2015, p.312). Appendix E specifies the methodology to generate an approximately Poisson-distributed random variable with finite range. Furthermore, we focus on order fulfilment systems with short lead times by assuming the lead time E of an order to be Bernoulli-distributed with parameter η in order to incorporate the customers' requirements for short order deadlines (see Chapter 1). Finally, we assume low variation for the processing performance per time unit L<sup>p</sup> of any process p ∈ P and model L<sup>p</sup> as Bernoulli-distributed random variable with parameter θ<sup>p</sup> and transformed range. Due to these assumptions, each stochastic parameter of the order fulfilment system is specified by one input variable in the Latin hypercube design.

## **C.2 Samples for Performance Analysis**

The Latin hypercube design for performance analysis of order fulfilment systems with pmax processes consists of (2 + 2 · pmax) input variables, one for each of the following parameters:


The ranges of the parameters are sample-specific. For the stochastic parameters, we differentiate between the range of the random variable and the range of the corresponding input variable of the Latin hypercube design. To ensure that the resulting order fulfilment systems correspond to reasonable system configurations, the chosen ranges are based on the following considerations: First, the maximum possible processing performance per time unit lp,max of any process p ∈ P corresponds to the minimum possible number of incoming orders per time period amin. Second, the minimum (maximum) possible value of the provided capacity cp,min (cp,max) at process p ∈ P is determined by the ratio of the minimum (maximum) possible number of incoming orders per time period amin (amax) to the maximum (minimum) possible processing performance per time unit lp,max (lp,min) of process p. The discretisation of the lead time E and the processing performances per time unit L<sup>p</sup> of the processes p ∈ P is constant in all samples since these parameters are Bernoulli-distributed. In contrast, the range of the number of incoming orders per time period A and thus its discretisation are sample-specific. Unstable order fulfilment systems (see equation (6.7)) are removed from the sample.

In this thesis, we generate six samples for performance analysis of order fulfilment systems, three of which consider single-stage systems (sample A1, A2, B) and two of which consider two-stage systems (sample C, D). Sample E contains three-stage systems. The specification of each sample is given in Table C.1.


Table C.1: Samples for performance analysis of order fulfilment systems.


## **C.3 Samples for Capacity Planning**

Samples for capacity planning differ from the ones for performance analysis regarding the input variables of the Latin hypercube design. For capacity planning, we do not require input variables for capacity cp, p ∈ P, but one input variable to specify the performance requirement. Thus, the Latin hypercube design for capacity planning in order fulfilment systems with pmax processes consists of (3 + pmax) input variables, one for each of the following parameters:


We determine the ranges of the parameters of the samples for capacity planning based on the same considerations as the ones of the samples for performance analysis (see Section C.2). To specify the performance requirement, any performance measure of the order fulfilment system (see Table 6.1) can be used. In this thesis, we focus on α-service level SLα. Thus, the range of the corresponding input variable of the Latin hypercube design is given by the interval [0,1].

We generate two samples (sample F, G) for capacity planning in order fulfilment systems, whereby sample F considers single-stage systems and sample G contains two-stage systems. The specification of each sample is given in Table C.2.


Table C.2: Samples for capacity planning in order fulfilment systems.

# **D Model Verification**

To verify the implementation of the exact model for performance analysis, we conduct a model comparison between the exact model and a simulation model. The simulation model depicts system behaviour of a multi-stage, stochastic order fulfilment system with levelled order release and customer-required deadlines, analogous to the Markov chain. Based on the simulation results, we calculate the same performance measures of the order fulfilment system as in the exact model (see Table 6.1). β-service level is used to determine the length of each simulation replication (precision of 1.0E-05) and the number of simulation replications (precision of 1.0E-04). Appendix B introduces the simulation model in detail.

As evaluation criteria, we calculate the absolute and relative deviations of selected performance measures – α-service level, β-service level, expected number of unprocessed orders, and expected number of completed orders – between the exact and the simulation model based on samples C and E. Sample C contains 256 two-stage order fulfilment systems, and sample E consists of 353 three-stage order fulfilment systems. The samples differ regarding the ranges of the parameters, especially regarding the discretisation of the number of incoming orders per time period A: Its range is given by A<sup>C</sup> = {3, 4, . . . , 14} in sample C and A<sup>E</sup> = {2, 3, . . . , 7} in sample E. A comprehensive description of the samples is given in Appendix C.2.

Table D.1 presents the absolute and relative deviations of the selected performance measures between the exact model and the simulation model in samples C and E. For each performance measure, Table D.1 gives the minimum, the maximum, and the average value, as well as the 2.5%- and the 97.5%-quantile of the absolute and the relative deviations, respectively.


Table D.1: Deviation of selected performance measures between the exact model and the simulation model.

<sup>1</sup> Difference between the value of the considered performance measure calculated based on the exact model and the one calculated based on the simulation model.

<sup>2</sup> Difference between the value of the considered performance measure calculated based on the exact and the one calculated based on the simulation model, divided by the value calculated based on the simulation model.

We observe absolute deviations of β-service level between -4.68E-04 and 3.96E-04 in sample C and between -1.67E-04 and 1.68E-04 in sample E. The corresponding 95%-confidence interval is given by [-6.98E-05, 3.16E-05] in sample C and [-6.99E-05, 5.89E-05] in sample E. Regarding α-service level, we observe absolute and relative deviations of the same magnitude as for β-service level. However, the absolute and relative deviations of α-service level are more fluctuating than those of β-service level as indicated by the wider 95%-confidence intervals of [-1.07E-04, 1.01.E-04] in sample C and [-1.37E-04, 1.24E-04] in sample E.

We observe absolute deviations of the expected number of unprocessed orders between -1.23E-02 and 3.10E-02 in sample C and between -4.61E-03 and 2.05E-02 in sample E. The corresponding 95%-confidence interval is given by [-6.24E-03, 2.02E-02] in sample C and [-2.5E-03, 1.03E-02] in sample E. The absolute and relative deviations of the expected number of completed orders are of the same magnitude as the ones of the expected number of unprocessed orders.

The deviations of the expected number of unprocessed orders and the expected number of completed orders are higher than the ones of α- and β-service level. However, we measure the expected number of unprocessed orders and the expected number of completed orders in number of orders, whereas α- and β-service level are calculated as a ratio with a range of [0,1]. Thus, for each performance measure, the deviations between the exact model and the simulation model are negligible despite the different magnitudes of the deviations. In conclusion, the results of the model comparison confirm that the exact model is implemented correctly.

# **E Generation of Approximately Poisson-distributed Random Variables with Finite Range**

Random variable K is said to be Poisson-distributed P oi(λ) with parameter λ > 0 and discrete range {0, 1, . . .} if its probability distribution is given by

$$P(K=k) = \begin{cases} \frac{1}{k!} \cdot e^{-\lambda} \cdot \lambda^k & k \in \{0, 1, \ldots\} \\ 0 & \text{otherwise.} \end{cases} \tag{\text{E.I.}}$$

The expected value and the variance are calculated as follows (Law 2015, p.312)

$$E(K) = \lambda \qquad \qquad Var(K) = \lambda. \tag{\mathbb{E}.2}$$

In the following, we present a methodology to generate the probability distribution of an approximately Poisson-distributed random variable J with finite range J = {jmin, . . . , jmax} whose expected value is given by E(J).

Initially, we calculate the parameter λ<sup>0</sup> of the corresponding Poisson-distributed random variable K as follows

$$
\lambda\_0 = E(K) = E(J) - j\_{min}.\tag{E.3}
$$

Subsequently, an iterative procedure starts. Each iteration i ∈ N consists of the following steps:

1. Determine the approximate upper bound Π<sup>i</sup> of the range of K ∼ P oi(λi−1) based on the following condition:

$$P(K \le \Pi\_i) \ge 1 - \epsilon\_P,\tag{E.4}$$

whereby <sup>P</sup> is set to 1.0E-09.

2. Calculate the scale parameter Φ<sup>i</sup> :

$$\Phi\_i = \frac{j\_{max} - j\_{min}}{\Pi\_i}.\tag{E.5}$$

3. Update the Poisson parameter λ<sup>i</sup> :

$$
\lambda\_i = \frac{E(J) - j\_{min}}{\Phi\_i}.\tag{E.6}
$$

4. Calculate the normalisation constant Ψ<sup>i</sup> based on K ∼ P oi(λi):

$$\Psi\_i = \sum\_{k=0}^{\Pi\_i} P(K=k). \tag{\mathbb{E}.7}$$

5. Check the stopping criterion:

$$
\Psi\_i > 1 - \epsilon\_P,\tag{E.8}
$$

whereby <sup>P</sup> is set to 1.0E-09. If the stopping criterion is met, the iterative procedure terminates with λ = λ<sup>i</sup> , Π = Π<sup>i</sup> , Φ = Φ<sup>i</sup> , and Ψ = Ψ<sup>i</sup> . Otherwise, iteration (i + 1) starts.

Finally, the probability distribution of J is derived from the probability distribution P oi(λ) of K as follows

$$\begin{aligned} P(J=j) &= \frac{1}{\Psi} \cdot \sum\_{k \in \mathbb{Z}(j)} P(K=k) & \forall j \in \mathcal{J} \\ \mathcal{Z}(j) &= \left\{ k \in \{0, \dots, \Pi\} \mid \left( j - \frac{1}{2} \right) \le \left( \Phi \cdot k + j\_{\min} \right) < \left( j + \frac{1}{2} \right) \right\} . \end{aligned} \tag{E.9}$$

# **Glossary of Notation**

### Description of an Order Fulfilment System





### Performance Measures of an Order Fulfilment System




### Discrete-time Markov Chain



### Analytical Model for Performance Analysis



### Derivative-free and Blackbox Optimisation



### Capacity Planning Problem


## Algorithmic Parameters of the MADS Algorithm


### Algorithmic Parameters of the SO-I Algorithm


### Simulation Model for Performance Analysis



### Design of Experiments



# **List of Figures**


# **List of Tables**





# **List of Publications**

Mohring, U., M. Baumann and K. Furmans (2020). Discrete-Time Analysis of Levelled Order Release and Staffing in Order Picking Systems. Logistics Research 13(9), p. 1.

# **References**


(eds.), Mathematics Without Boundaries: Surveys in Interdisciplinary Research, p. 31–56. New York, NY: Springer New York.


Yu, M. and R. B. De Koster (2009). The impact of order batching and picking area zoning on order picking system performance. European Journal of Operational Research 198(2), p. 480–490.

Wissenschaftliche Berichte des Instituts für Fördertechnik und Logistiksysteme des Karlsruher Instituts für Technologie (KIT) Prof. Dr.-Ing. Kai Furmans [Hrsg.]

Order fulfilment systems are forced to efficiently manage a highly volatile customer demand while simultaneously meeting customer-required short order deadlines. To handle these challenges, we develop the *Strategy of Levelled Order Release* (*LOR*) that balances workload in multi-stage stochastic systems with order deadlines over time. The contribution of this work is three-fold: We introduce (1) the workload balancing concept *LOR*, (2) a discrete-time Markov chain for performance analysis, and (3) an algorithm for capacity planning under performance constraints in order fulfilment systems with *LOR*.

The Markov chain models order fulfilment according to *LOR* in multi-stage stochastic systems with customer-required order deadlines. It enables the exact calculation of multiple key figures, e.g. system throughput, service level, and backlog duration, for systems with general distributed input parameters. The capacity planning problem determines the minimum capacity that is required to guarantee specific performance requirements of the customers. As it is a blackbox optimisation problem, we solve it using the derivative-free optimisation algorithms *Mesh Adaptive Direct Search* and *Surrogate Optimisation Integer*. Numerical studies show that *LOR* achieves higher service levels and reduces the amount of required capacity compared to *FCFS*, especially in systems with high utilisation.

ISSN 0171-2772 ISBN 978-3-7315-1211-0